From 65401b57e5fedbce4d2b2d7bae832f4ec5c99872 Mon Sep 17 00:00:00 2001 From: Timo Betcke Date: Sun, 3 Nov 2024 18:36:58 +0000 Subject: [PATCH] Refactored boundary assembler --- benches/assembly_benchmark.rs | 26 +- examples/assembly.rs | 97 +- examples/test_parallel_assembly.rs | 509 ++-- src/assembly.rs | 453 ++-- src/assembly/boundary.rs | 31 +- src/assembly/boundary/assemblers.rs | 46 +- .../assemblers/adjoint_double_layer.rs | 46 - .../boundary/assemblers/double_layer.rs | 38 - .../boundary/assemblers/hypersingular.rs | 60 - .../boundary/assemblers/single_layer.rs | 36 - src/bindings.rs | 2138 ++++++++--------- src/helmholtz.rs | 146 ++ src/laplace.rs | 126 + src/lib.rs | 2 + tests/compare_to_bempp_cl.rs | 81 +- 15 files changed, 1972 insertions(+), 1863 deletions(-) delete mode 100644 src/assembly/boundary/assemblers/adjoint_double_layer.rs delete mode 100644 src/assembly/boundary/assemblers/double_layer.rs delete mode 100644 src/assembly/boundary/assemblers/hypersingular.rs delete mode 100644 src/assembly/boundary/assemblers/single_layer.rs create mode 100644 src/helmholtz.rs create mode 100644 src/laplace.rs diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs index 815de35b..df162635 100644 --- a/benches/assembly_benchmark.rs +++ b/benches/assembly_benchmark.rs @@ -1,7 +1,11 @@ -use bempp::assembly::boundary::BoundaryAssembler; +use bempp::assembly::boundary::integrands::SingleLayerBoundaryIntegrand; +use bempp::assembly::boundary::{BoundaryAssembler, BoundaryAssemblerOptions}; +use bempp::assembly::kernels::KernelEvaluator; use bempp::function::SerialFunctionSpace; use bempp::traits::{BoundaryAssembly, FunctionSpace}; use criterion::{criterion_group, criterion_main, Criterion}; +use green_kernels::laplace_3d::Laplace3dKernel; +use green_kernels::types::GreenKernelEvalType; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::{Continuity, ReferenceCellType}; use ndgrid::shapes::regular_sphere; @@ -19,13 +23,21 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) { let mut matrix = rlst_dynamic_array2!(f64, [space.global_size(), space.global_size()]); let colouring = space.cell_colouring(); - let mut a = BoundaryAssembler::::new_laplace_single_layer(); - a.set_quadrature_degree(ReferenceCellType::Triangle, 16); - a.set_singular_quadrature_degree( + + let mut assembler = BoundaryAssembler::::new( + SingleLayerBoundaryIntegrand::new(), + KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value), + BoundaryAssemblerOptions::default(), + 1, + 0, + ); + + assembler.set_quadrature_degree(ReferenceCellType::Triangle, 16); + assembler.set_singular_quadrature_degree( (ReferenceCellType::Triangle, ReferenceCellType::Triangle), 4, ); - a.set_batch_size(128); + assembler.set_batch_size(128); group.bench_function( format!( @@ -33,7 +45,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) { space.global_size(), space.global_size() ), - |b| b.iter(|| a.assemble_singular_into_dense(&mut matrix, &space, &space)), + |b| b.iter(|| assembler.assemble_singular_into_dense(&mut matrix, &space, &space)), ); group.bench_function( format!( @@ -43,7 +55,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) { ), |b| { b.iter(|| { - a.assemble_nonsingular_into_dense( + assembler.assemble_nonsingular_into_dense( &mut matrix, &space, &space, diff --git a/examples/assembly.rs b/examples/assembly.rs index 0f94e8d4..4454ecb1 100644 --- a/examples/assembly.rs +++ b/examples/assembly.rs @@ -1,10 +1,11 @@ -use bempp::assembly::boundary::BoundaryAssembler; +use bempp::assembly::boundary::{BoundaryAssembler, BoundaryAssemblerOptions}; use bempp::function::SerialFunctionSpace; +use bempp::laplace::assembler::laplace_single_layer; use bempp::traits::{BoundaryAssembly, FunctionSpace}; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::{Continuity, ReferenceCellType}; use ndgrid::shapes::regular_sphere; -use rlst::{rlst_dynamic_array2, RandomAccessByRef}; +use rlst::{rlst_dynamic_array2, RandomAccessByRef, Shape}; fn main() { // Create a grid, family of elements, and function space @@ -12,73 +13,67 @@ fn main() { let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = SerialFunctionSpace::new(&grid, &element); - // Create an array to store the assembled discrete operator - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - - // Create an assembler for the Laplace single layer operator - let mut a = BoundaryAssembler::::new_laplace_single_layer(); - // Adjust the quadrature degree for non-singular integrals on a triangle. // This makes the integrals use a quadrature rule with 16 points - a.set_quadrature_degree(ReferenceCellType::Triangle, 16); + let mut options = BoundaryAssemblerOptions::default(); + options.set_regular_quadrature_degree(ReferenceCellType::Triangle, 16); // Adjust the quadrature degree for singular integrals on a pair ortriangles. // This makes the integrals use a quadrature rule based on a rule on an interval with 4 points - a.set_singular_quadrature_degree( + options.set_singular_quadrature_degree( (ReferenceCellType::Triangle, ReferenceCellType::Triangle), 4, ); - // Assemble the discrete operator - a.assemble_into_dense(&mut matrix, &space, &space); + // Assemble the single layer Laplace boundary operator. + let matrix = laplace_single_layer(&space, &space, options.clone()); // Print the entries of the matrix println!("Lagrange single layer matrix"); - for i in 0..ndofs { - for j in 0..ndofs { + for i in 0..matrix.shape()[0] { + for j in 0..matrix.shape()[1] { print!("{:.4} ", matrix.get([i, j]).unwrap()); } println!(); } println!(); - // Assemble just the singular part - let mut singular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - a.assemble_singular_into_dense(&mut singular_matrix, &space, &space); - println!("Lagrange single layer matrix (singular part)"); - for i in 0..ndofs { - for j in 0..ndofs { - print!("{:.4} ", singular_matrix.get([i, j]).unwrap()); - } - println!(); - } - println!(); + // // Assemble just the singular part + // let mut singular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); + // a.assemble_singular_into_dense(&mut singular_matrix, &space, &space); + // println!("Lagrange single layer matrix (singular part)"); + // for i in 0..ndofs { + // for j in 0..ndofs { + // print!("{:.4} ", singular_matrix.get([i, j]).unwrap()); + // } + // println!(); + // } + // println!(); - // For grids with a larger number of cells, the singular part will be sparse. It can be assembled as a CSR matrix as follows - println!("Lagrange single layer matrix (singular part) as CSR matrix"); - let singular_sparse_matrix = a.assemble_singular_into_csr(&space, &space); - println!("indices: {:?}", singular_sparse_matrix.indices()); - println!("indptr: {:?}", singular_sparse_matrix.indptr()); - println!("data: {:?}", singular_sparse_matrix.data()); - println!(); + // // For grids with a larger number of cells, the singular part will be sparse. It can be assembled as a CSR matrix as follows + // println!("Lagrange single layer matrix (singular part) as CSR matrix"); + // let singular_sparse_matrix = a.assemble_singular_into_csr(&space, &space); + // println!("indices: {:?}", singular_sparse_matrix.indices()); + // println!("indptr: {:?}", singular_sparse_matrix.indptr()); + // println!("data: {:?}", singular_sparse_matrix.data()); + // println!(); - // Assemble just the non-singular part - let colouring = space.cell_colouring(); - let mut nonsingular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - a.assemble_nonsingular_into_dense( - &mut nonsingular_matrix, - &space, - &space, - &colouring, - &colouring, - ); - println!("Lagrange single layer matrix (nonsingular part)"); - for i in 0..ndofs { - for j in 0..ndofs { - print!("{:.4} ", nonsingular_matrix.get([i, j]).unwrap()); - } - println!(); - } - println!(); + // // Assemble just the non-singular part + // let colouring = space.cell_colouring(); + // let mut nonsingular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); + // a.assemble_nonsingular_into_dense( + // &mut nonsingular_matrix, + // &space, + // &space, + // &colouring, + // &colouring, + // ); + // println!("Lagrange single layer matrix (nonsingular part)"); + // for i in 0..ndofs { + // for j in 0..ndofs { + // print!("{:.4} ", nonsingular_matrix.get([i, j]).unwrap()); + // } + // println!(); + // } + // println!(); } diff --git a/examples/test_parallel_assembly.rs b/examples/test_parallel_assembly.rs index bb3b2f1a..bb6d3bcb 100644 --- a/examples/test_parallel_assembly.rs +++ b/examples/test_parallel_assembly.rs @@ -1,255 +1,256 @@ -//? mpirun -n {{NPROCESSES}} --features "mpi" - -#[cfg(feature = "mpi")] -use approx::assert_relative_eq; -#[cfg(feature = "mpi")] -use bempp::{ - assembly::boundary::BoundaryAssembler, - function::{ParallelFunctionSpace, SerialFunctionSpace}, - traits::{BoundaryAssembly, FunctionSpace, ParallelBoundaryAssembly}, -}; -#[cfg(feature = "mpi")] -use itertools::izip; -#[cfg(feature = "mpi")] -use mpi::{ - collective::CommunicatorCollectives, - environment::Universe, - request::WaitGuard, - traits::{Communicator, Destination, Source}, -}; -#[cfg(feature = "mpi")] -use ndelement::{ - ciarlet::{CiarletElement, LagrangeElementFamily}, - types::{Continuity, ReferenceCellType}, -}; -#[cfg(feature = "mpi")] -use ndgrid::{ - grid::parallel::ParallelGrid, - traits::{Builder, Entity, Grid, ParallelBuilder}, - SingleElementGrid, SingleElementGridBuilder, -}; -#[cfg(feature = "mpi")] -use rlst::{CsrMatrix, Shape}; -#[cfg(feature = "mpi")] -use std::collections::{hash_map::Entry, HashMap}; - -#[cfg(feature = "mpi")] -fn create_single_element_grid_data(b: &mut SingleElementGridBuilder, n: usize) { - for y in 0..n { - for x in 0..n { - b.add_point( - y * n + x, - &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64, 0.0], - ); - } - } - - for i in 0..n - 1 { - for j in 0..n - 1 { - b.add_cell( - i * (n - 1) + j, - &[j * n + i, j * n + i + 1, j * n + i + n, j * n + i + n + 1], - ); - } - } -} - -#[cfg(feature = "mpi")] -fn example_single_element_grid( - comm: &C, - n: usize, -) -> ParallelGrid<'_, C, SingleElementGrid>> { - let rank = comm.rank(); - - let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); - - if rank == 0 { - create_single_element_grid_data(&mut b, n); - b.create_parallel_grid(comm) - } else { - b.receive_parallel_grid(comm, 0) - } -} - -#[cfg(feature = "mpi")] -fn example_single_element_grid_serial(n: usize) -> SingleElementGrid> { - let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); - create_single_element_grid_data(&mut b, n); - b.create_grid() -} - -#[cfg(feature = "mpi")] -fn test_parallel_assembly_single_element_grid( - comm: &C, - degree: usize, - cont: Continuity, -) { - let rank = comm.rank(); - let size = comm.size(); - - let n = 10; - let grid = example_single_element_grid(comm, n); - let element = LagrangeElementFamily::::new(degree, cont); - let space = ParallelFunctionSpace::new(&grid, &element); - - let a = BoundaryAssembler::::new_laplace_single_layer(); - - let matrix = a.parallel_assemble_singular_into_csr(&space, &space); - - if rank == 0 { - // Compute the same matrix on a single process - let serial_grid = example_single_element_grid_serial(n); - let serial_space = SerialFunctionSpace::new(&serial_grid, &element); - let serial_matrix = a.assemble_singular_into_csr(&serial_space, &serial_space); - - // Dofs associated with each cell (by cell id) - let mut serial_dofmap = HashMap::new(); - for cell in serial_grid.entity_iter(2) { - serial_dofmap.insert( - cell.id().unwrap(), - serial_space - .cell_dofs(cell.local_index()) - .unwrap() - .iter() - .map(|i| serial_space.global_dof_index(*i)) - .collect::>(), - ); - } - let mut parallel_dofmap = HashMap::new(); - for cell in grid.entity_iter(2) { - parallel_dofmap.insert( - cell.id().unwrap(), - space - .cell_dofs(cell.local_index()) - .unwrap() - .iter() - .map(|i| space.global_dof_index(*i)) - .collect::>(), - ); - } - for p in 1..size { - let process = comm.process_at_rank(p); - let (cell_ids, _status) = process.receive_vec::(); - let (dofs, _status) = process.receive_vec::(); - let (dofs_len, _status) = process.receive_vec::(); - let mut start = 0; - for (id, len) in izip!(cell_ids, dofs_len) { - if let Entry::Vacant(e) = parallel_dofmap.entry(id) { - e.insert(dofs[start..start + len].to_vec()); - } else { - assert_eq!(parallel_dofmap[&id], dofs[start..start + len]); - } - start += len; - } - } - - let mut index_map = vec![0; serial_space.global_size()]; - - for (id, dofs) in parallel_dofmap { - for (i, j) in izip!(&serial_dofmap[&id], dofs) { - index_map[j] = *i; - } - } - - // Gather sparse matrices onto process 0 - let mut rows = vec![]; - let mut cols = vec![]; - let mut data = vec![]; - - let mut r = 0; - for (i, index) in matrix.indices().iter().enumerate() { - while i >= matrix.indptr()[r + 1] { - r += 1; - } - rows.push(index_map[r]); - cols.push(index_map[*index]); - data.push(matrix.data()[i]); - } - - for p in 1..size { - let process = comm.process_at_rank(p); - let (indices, _status) = process.receive_vec::(); - let (indptr, _status) = process.receive_vec::(); - let (subdata, _status) = process.receive_vec::(); - let mat = CsrMatrix::new(matrix.shape(), indices, indptr, subdata); - - let mut r = 0; - for (i, index) in mat.indices().iter().enumerate() { - while i >= mat.indptr()[r + 1] { - r += 1; - } - rows.push(index_map[r]); - cols.push(index_map[*index]); - data.push(mat.data()[i]); - } - } - let full_matrix = CsrMatrix::from_aij( - [space.global_size(), space.global_size()], - &rows, - &cols, - &data, - ) - .unwrap(); - - // Compare to matrix assembled on just this process - for (i, j) in full_matrix.indices().iter().zip(serial_matrix.indices()) { - assert_eq!(i, j); - } - for (i, j) in full_matrix.indptr().iter().zip(serial_matrix.indptr()) { - assert_eq!(i, j); - } - for (i, j) in full_matrix.data().iter().zip(serial_matrix.data()) { - assert_relative_eq!(i, j, epsilon = 1e-10); - } - } else { - let mut cell_ids = vec![]; - let mut dofs = vec![]; - let mut dofs_len = vec![]; - for cell in grid.entity_iter(2) { - cell_ids.push(cell.id().unwrap()); - let cell_dofs = space - .cell_dofs(cell.local_index()) - .unwrap() - .iter() - .map(|i| space.global_dof_index(*i)) - .collect::>(); - dofs.extend_from_slice(&cell_dofs); - dofs_len.push(cell_dofs.len()); - } - - mpi::request::scope(|scope| { - let root = comm.process_at_rank(0); - // TODO: send this: - let _ = WaitGuard::from(root.immediate_send(scope, &cell_ids)); - let _ = WaitGuard::from(root.immediate_send(scope, &dofs)); - let _ = WaitGuard::from(root.immediate_send(scope, &dofs_len)); - let _ = WaitGuard::from(root.immediate_send(scope, matrix.indices())); - let _ = WaitGuard::from(root.immediate_send(scope, matrix.indptr())); - let _ = WaitGuard::from(root.process_at_rank(0).immediate_send(scope, matrix.data())); - }); - } -} - -#[cfg(feature = "mpi")] -fn main() { - let universe: Universe = mpi::initialize().unwrap(); - let world = universe.world(); - let rank = world.rank(); - - for degree in 0..4 { - if rank == 0 { - println!("Testing assembly with DP{degree} using SingleElementGrid in parallel."); - } - test_parallel_assembly_single_element_grid(&world, degree, Continuity::Discontinuous); - world.barrier(); - } - for degree in 1..4 { - if rank == 0 { - println!("Testing assembly with P{degree} using SingleElementGrid in parallel."); - } - test_parallel_assembly_single_element_grid(&world, degree, Continuity::Standard); - world.barrier(); - } -} -#[cfg(not(feature = "mpi"))] +// //? mpirun -n {{NPROCESSES}} --features "mpi" + +// #[cfg(feature = "mpi")] +// use approx::assert_relative_eq; +// #[cfg(feature = "mpi")] +// use bempp::{ +// assembly::boundary::BoundaryAssembler, +// function::{ParallelFunctionSpace, SerialFunctionSpace}, +// traits::{BoundaryAssembly, FunctionSpace, ParallelBoundaryAssembly}, +// }; +// #[cfg(feature = "mpi")] +// use itertools::izip; +// #[cfg(feature = "mpi")] +// use mpi::{ +// collective::CommunicatorCollectives, +// environment::Universe, +// request::WaitGuard, +// traits::{Communicator, Destination, Source}, +// }; +// #[cfg(feature = "mpi")] +// use ndelement::{ +// ciarlet::{CiarletElement, LagrangeElementFamily}, +// types::{Continuity, ReferenceCellType}, +// }; +// #[cfg(feature = "mpi")] +// use ndgrid::{ +// grid::parallel::ParallelGrid, +// traits::{Builder, Entity, Grid, ParallelBuilder}, +// SingleElementGrid, SingleElementGridBuilder, +// }; +// #[cfg(feature = "mpi")] +// use rlst::{CsrMatrix, Shape}; +// #[cfg(feature = "mpi")] +// use std::collections::{hash_map::Entry, HashMap}; + +// #[cfg(feature = "mpi")] +// fn create_single_element_grid_data(b: &mut SingleElementGridBuilder, n: usize) { +// for y in 0..n { +// for x in 0..n { +// b.add_point( +// y * n + x, +// &[x as f64 / (n - 1) as f64, y as f64 / (n - 1) as f64, 0.0], +// ); +// } +// } + +// for i in 0..n - 1 { +// for j in 0..n - 1 { +// b.add_cell( +// i * (n - 1) + j, +// &[j * n + i, j * n + i + 1, j * n + i + n, j * n + i + n + 1], +// ); +// } +// } +// } + +// #[cfg(feature = "mpi")] +// fn example_single_element_grid( +// comm: &C, +// n: usize, +// ) -> ParallelGrid<'_, C, SingleElementGrid>> { +// let rank = comm.rank(); + +// let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); + +// if rank == 0 { +// create_single_element_grid_data(&mut b, n); +// b.create_parallel_grid(comm) +// } else { +// b.receive_parallel_grid(comm, 0) +// } +// } + +// #[cfg(feature = "mpi")] +// fn example_single_element_grid_serial(n: usize) -> SingleElementGrid> { +// let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); +// create_single_element_grid_data(&mut b, n); +// b.create_grid() +// } + +// // #[cfg(feature = "mpi")] +// // fn test_parallel_assembly_single_element_grid( +// // comm: &C, +// // degree: usize, +// // cont: Continuity, +// // ) { +// // let rank = comm.rank(); +// // let size = comm.size(); + +// // let n = 10; +// // let grid = example_single_element_grid(comm, n); +// // let element = LagrangeElementFamily::::new(degree, cont); +// // let space = ParallelFunctionSpace::new(&grid, &element); + +// // let a = BoundaryAssembler::::new_laplace_single_layer(); + +// // let matrix = a.parallel_assemble_singular_into_csr(&space, &space); + +// // if rank == 0 { +// // // Compute the same matrix on a single process +// // let serial_grid = example_single_element_grid_serial(n); +// // let serial_space = SerialFunctionSpace::new(&serial_grid, &element); +// // let serial_matrix = a.assemble_singular_into_csr(&serial_space, &serial_space); + +// // // Dofs associated with each cell (by cell id) +// // let mut serial_dofmap = HashMap::new(); +// // for cell in serial_grid.entity_iter(2) { +// // serial_dofmap.insert( +// // cell.id().unwrap(), +// // serial_space +// // .cell_dofs(cell.local_index()) +// // .unwrap() +// // .iter() +// // .map(|i| serial_space.global_dof_index(*i)) +// // .collect::>(), +// // ); +// // } +// // let mut parallel_dofmap = HashMap::new(); +// // for cell in grid.entity_iter(2) { +// // parallel_dofmap.insert( +// // cell.id().unwrap(), +// // space +// // .cell_dofs(cell.local_index()) +// // .unwrap() +// // .iter() +// // .map(|i| space.global_dof_index(*i)) +// // .collect::>(), +// // ); +// // } +// // for p in 1..size { +// // let process = comm.process_at_rank(p); +// // let (cell_ids, _status) = process.receive_vec::(); +// // let (dofs, _status) = process.receive_vec::(); +// // let (dofs_len, _status) = process.receive_vec::(); +// // let mut start = 0; +// // for (id, len) in izip!(cell_ids, dofs_len) { +// // if let Entry::Vacant(e) = parallel_dofmap.entry(id) { +// // e.insert(dofs[start..start + len].to_vec()); +// // } else { +// // assert_eq!(parallel_dofmap[&id], dofs[start..start + len]); +// // } +// // start += len; +// // } +// // } + +// // let mut index_map = vec![0; serial_space.global_size()]; + +// // for (id, dofs) in parallel_dofmap { +// // for (i, j) in izip!(&serial_dofmap[&id], dofs) { +// // index_map[j] = *i; +// // } +// // } + +// // // Gather sparse matrices onto process 0 +// // let mut rows = vec![]; +// // let mut cols = vec![]; +// // let mut data = vec![]; + +// // let mut r = 0; +// // for (i, index) in matrix.indices().iter().enumerate() { +// // while i >= matrix.indptr()[r + 1] { +// // r += 1; +// // } +// // rows.push(index_map[r]); +// // cols.push(index_map[*index]); +// // data.push(matrix.data()[i]); +// // } + +// // for p in 1..size { +// // let process = comm.process_at_rank(p); +// // let (indices, _status) = process.receive_vec::(); +// // let (indptr, _status) = process.receive_vec::(); +// // let (subdata, _status) = process.receive_vec::(); +// // let mat = CsrMatrix::new(matrix.shape(), indices, indptr, subdata); + +// // let mut r = 0; +// // for (i, index) in mat.indices().iter().enumerate() { +// // while i >= mat.indptr()[r + 1] { +// // r += 1; +// // } +// // rows.push(index_map[r]); +// // cols.push(index_map[*index]); +// // data.push(mat.data()[i]); +// // } +// // } +// // let full_matrix = CsrMatrix::from_aij( +// // [space.global_size(), space.global_size()], +// // &rows, +// // &cols, +// // &data, +// // ) +// // .unwrap(); + +// // // Compare to matrix assembled on just this process +// // for (i, j) in full_matrix.indices().iter().zip(serial_matrix.indices()) { +// // assert_eq!(i, j); +// // } +// // for (i, j) in full_matrix.indptr().iter().zip(serial_matrix.indptr()) { +// // assert_eq!(i, j); +// // } +// // for (i, j) in full_matrix.data().iter().zip(serial_matrix.data()) { +// // assert_relative_eq!(i, j, epsilon = 1e-10); +// // } +// // } else { +// // let mut cell_ids = vec![]; +// // let mut dofs = vec![]; +// // let mut dofs_len = vec![]; +// // for cell in grid.entity_iter(2) { +// // cell_ids.push(cell.id().unwrap()); +// // let cell_dofs = space +// // .cell_dofs(cell.local_index()) +// // .unwrap() +// // .iter() +// // .map(|i| space.global_dof_index(*i)) +// // .collect::>(); +// // dofs.extend_from_slice(&cell_dofs); +// // dofs_len.push(cell_dofs.len()); +// // } + +// // mpi::request::scope(|scope| { +// // let root = comm.process_at_rank(0); +// // // TODO: send this: +// // let _ = WaitGuard::from(root.immediate_send(scope, &cell_ids)); +// // let _ = WaitGuard::from(root.immediate_send(scope, &dofs)); +// // let _ = WaitGuard::from(root.immediate_send(scope, &dofs_len)); +// // let _ = WaitGuard::from(root.immediate_send(scope, matrix.indices())); +// // let _ = WaitGuard::from(root.immediate_send(scope, matrix.indptr())); +// // let _ = WaitGuard::from(root.process_at_rank(0).immediate_send(scope, matrix.data())); +// // }); +// // } +// // } + +// // #[cfg(feature = "mpi")] +// // fn main() { +// // let universe: Universe = mpi::initialize().unwrap(); +// // let world = universe.world(); +// // let rank = world.rank(); + +// // for degree in 0..4 { +// // if rank == 0 { +// // println!("Testing assembly with DP{degree} using SingleElementGrid in parallel."); +// // } +// // test_parallel_assembly_single_element_grid(&world, degree, Continuity::Discontinuous); +// // world.barrier(); +// // } +// // for degree in 1..4 { +// // if rank == 0 { +// // println!("Testing assembly with P{degree} using SingleElementGrid in parallel."); +// // } +// // test_parallel_assembly_single_element_grid(&world, degree, Continuity::Standard); +// // world.barrier(); +// // } +// // } +// fn main() {} + fn main() {} diff --git a/src/assembly.rs b/src/assembly.rs index 9e39be3a..ba2ad65d 100644 --- a/src/assembly.rs +++ b/src/assembly.rs @@ -5,240 +5,239 @@ pub mod fmm_tools; pub mod kernels; pub mod potential; -#[cfg(test)] -mod test { - use super::*; - use crate::function::SerialFunctionSpace; - use crate::traits::{BoundaryAssembly, FunctionSpace}; - use cauchy::{c32, c64}; - use ndelement::ciarlet::CiarletElement; - use ndelement::ciarlet::LagrangeElementFamily; - use ndelement::types::{Continuity, ReferenceCellType}; - use ndgrid::{ - grid::serial::{SingleElementGrid, SingleElementGridBuilder}, - shapes::regular_sphere, - traits::Builder, - types::RealScalar, - }; - use paste::paste; - use rlst::{rlst_dynamic_array2, MatrixInverse, RlstScalar}; +// #[cfg(test)] +// mod test { +// use super::*; +// use crate::function::SerialFunctionSpace; +// use crate::traits::{BoundaryAssembly, FunctionSpace}; +// use cauchy::{c32, c64}; +// use ndelement::ciarlet::CiarletElement; +// use ndelement::ciarlet::LagrangeElementFamily; +// use ndelement::types::{Continuity, ReferenceCellType}; +// use ndgrid::{ +// grid::serial::{SingleElementGrid, SingleElementGridBuilder}, +// shapes::regular_sphere, +// traits::Builder, +// types::RealScalar, +// }; +// use paste::paste; +// use rlst::{rlst_dynamic_array2, MatrixInverse, RlstScalar}; - fn quadrilateral_grid() -> SingleElementGrid> - { - let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); - for j in 0..4 { - for i in 0..4 { - b.add_point( - 4 * j + i, - &[ - num::cast::(i).unwrap() / num::cast::(3.0).unwrap(), - num::cast::(j).unwrap() / num::cast::(3.0).unwrap(), - num::cast::(0.0).unwrap(), - ], - ); - } - } - for j in 0..3 { - for i in 0..3 { - b.add_cell( - 3 * j + i, - &[4 * j + i, 4 * j + i + 1, 4 * j + i + 4, 4 * j + i + 5], - ); - } - } - b.create_grid() - } +// fn quadrilateral_grid() -> SingleElementGrid> +// { +// let mut b = SingleElementGridBuilder::::new(3, (ReferenceCellType::Quadrilateral, 1)); +// for j in 0..4 { +// for i in 0..4 { +// b.add_point( +// 4 * j + i, +// &[ +// num::cast::(i).unwrap() / num::cast::(3.0).unwrap(), +// num::cast::(j).unwrap() / num::cast::(3.0).unwrap(), +// num::cast::(0.0).unwrap(), +// ], +// ); +// } +// } +// for j in 0..3 { +// for i in 0..3 { +// b.add_cell( +// 3 * j + i, +// &[4 * j + i, 4 * j + i + 1, 4 * j + i + 4, 4 * j + i + 5], +// ); +// } +// } +// b.create_grid() +// } - /* - fn mixed_grid>() -> MixedGrid - where - for<'a> Array, 2>, 2>, 2>: - MatrixInverse, - { - let mut b = MixedGridBuilder::<3, T>::new(()); - for j in 0..4 { - for i in 0..4 { - b.add_point( - 4 * j + i, - [ - num::cast::(i).unwrap() / num::cast::(3.0).unwrap(), - num::cast::(j).unwrap() / num::cast::(3.0).unwrap(), - num::cast::(0.0).unwrap(), - ], - ); - } - } - for j in 0..3 { - b.add_cell( - j, - ( - vec![4 * j, 4 * j + 1, 4 * j + 4, 4 * j + 5], - ReferenceCellType::Quadrilateral, - 1, - ), - ); - } - for j in 0..3 { - b.add_cell( - 3 + 2 * j, - ( - vec![4 * j + 1, 4 * j + 2, 4 * j + 6], - ReferenceCellType::Triangle, - 1, - ), - ); - b.add_cell( - 4 + 2 * j, - ( - vec![4 * j + 1, 4 * j + 6, 4 * j + 5], - ReferenceCellType::Triangle, - 1, - ), - ); - } - for j in 0..3 { - b.add_cell( - 9 + j, - ( - vec![4 * j + 2, 4 * j + 3, 4 * j + 6, 4 * j + 7], - ReferenceCellType::Quadrilateral, - 1, - ), - ); - } - b.create_grid() - } - */ +// /* +// fn mixed_grid>() -> MixedGrid +// where +// for<'a> Array, 2>, 2>, 2>: +// MatrixInverse, +// { +// let mut b = MixedGridBuilder::<3, T>::new(()); +// for j in 0..4 { +// for i in 0..4 { +// b.add_point( +// 4 * j + i, +// [ +// num::cast::(i).unwrap() / num::cast::(3.0).unwrap(), +// num::cast::(j).unwrap() / num::cast::(3.0).unwrap(), +// num::cast::(0.0).unwrap(), +// ], +// ); +// } +// } +// for j in 0..3 { +// b.add_cell( +// j, +// ( +// vec![4 * j, 4 * j + 1, 4 * j + 4, 4 * j + 5], +// ReferenceCellType::Quadrilateral, +// 1, +// ), +// ); +// } +// for j in 0..3 { +// b.add_cell( +// 3 + 2 * j, +// ( +// vec![4 * j + 1, 4 * j + 2, 4 * j + 6], +// ReferenceCellType::Triangle, +// 1, +// ), +// ); +// b.add_cell( +// 4 + 2 * j, +// ( +// vec![4 * j + 1, 4 * j + 6, 4 * j + 5], +// ReferenceCellType::Triangle, +// 1, +// ), +// ); +// } +// for j in 0..3 { +// b.add_cell( +// 9 + j, +// ( +// vec![4 * j + 2, 4 * j + 3, 4 * j + 6, 4 * j + 7], +// ReferenceCellType::Quadrilateral, +// 1, +// ), +// ); +// } +// b.create_grid() +// } +// */ +// macro_rules! example_grid { +// (Triangle, $dtype:ident) => { +// regular_sphere(0) +// }; +// (Quadrilateral, $dtype:ident) => { +// quadrilateral_grid::<<$dtype as RlstScalar>::Real>() +// }; //(Mixed, $dtype:ident) => { +// // mixed_grid::<<$dtype as RlstScalar>::Real>() +// //}; +// } +// macro_rules! create_assembler { +// (Laplace, Hypersingular, $dtype:ident) => { +// paste! { +// boundary::HypersingularAssembler::<[<$dtype>], _, _>::new_laplace() +// } +// }; +// (Helmholtz, Hypersingular, $dtype:ident) => { +// paste! { +// boundary::HypersingularAssembler::<[<$dtype>], _, _>::new_helmholtz(3.0) +// } +// }; +// (Laplace, $operator:ident, $dtype:ident) => { +// paste! { +// boundary::BoundaryAssembler::<[<$dtype>], _, _>::[]() +// } +// }; +// (Helmholtz, $operator:ident, $dtype:ident) => { +// paste! { +// boundary::BoundaryAssembler::<[<$dtype>], _, _>::[](3.0) +// } +// }; +// } +// macro_rules! test_assembly { - macro_rules! example_grid { - (Triangle, $dtype:ident) => { - regular_sphere(0) - }; - (Quadrilateral, $dtype:ident) => { - quadrilateral_grid::<<$dtype as RlstScalar>::Real>() - }; //(Mixed, $dtype:ident) => { - // mixed_grid::<<$dtype as RlstScalar>::Real>() - //}; - } - macro_rules! create_assembler { - (Laplace, Hypersingular, $dtype:ident) => { - paste! { - boundary::HypersingularAssembler::<[<$dtype>], _, _>::new_laplace() - } - }; - (Helmholtz, Hypersingular, $dtype:ident) => { - paste! { - boundary::HypersingularAssembler::<[<$dtype>], _, _>::new_helmholtz(3.0) - } - }; - (Laplace, $operator:ident, $dtype:ident) => { - paste! { - boundary::BoundaryAssembler::<[<$dtype>], _, _>::[]() - } - }; - (Helmholtz, $operator:ident, $dtype:ident) => { - paste! { - boundary::BoundaryAssembler::<[<$dtype>], _, _>::[](3.0) - } - }; - } - macro_rules! test_assembly { +// ($(($dtype:ident, $pde:ident, $operator:ident, $cell:ident)),+) => { - ($(($dtype:ident, $pde:ident, $operator:ident, $cell:ident)),+) => { +// $( +// paste! { - $( - paste! { +// #[test] +// fn []() { - #[test] - fn []() { +// let grid = example_grid!($cell, $dtype); +// let element = LagrangeElementFamily::<[<$dtype>]>::new(0, Continuity::Discontinuous); +// let space = SerialFunctionSpace::new(&grid, &element); - let grid = example_grid!($cell, $dtype); - let element = LagrangeElementFamily::<[<$dtype>]>::new(0, Continuity::Discontinuous); - let space = SerialFunctionSpace::new(&grid, &element); +// let ndofs = space.global_size(); +// let mut matrix = rlst_dynamic_array2!([<$dtype>], [ndofs, ndofs]); - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!([<$dtype>], [ndofs, ndofs]); +// let a = create_assembler!($pde, $operator, $dtype); +// a.assemble_into_dense(&mut matrix, &space, &space); +// } - let a = create_assembler!($pde, $operator, $dtype); - a.assemble_into_dense(&mut matrix, &space, &space); - } +// } +// )* +// }; +// } - } - )* - }; - } - - test_assembly!( - (f64, Laplace, single_layer, Triangle), - (f32, Laplace, single_layer, Triangle), - //(c64, Laplace, single_layer, Triangle), - //(c32, Laplace, single_layer, Triangle), - (f64, Laplace, double_layer, Triangle), - (f32, Laplace, double_layer, Triangle), - //(c64, Laplace, double_layer, Triangle), - //(c32, Laplace, double_layer, Triangle), - (f64, Laplace, adjoint_double_layer, Triangle), - (f32, Laplace, adjoint_double_layer, Triangle), - //(c64, Laplace, adjoint_double_layer, Triangle), - //(c32, Laplace, adjoint_double_layer, Triangle), - (f64, Laplace, hypersingular, Triangle), - (f32, Laplace, hypersingular, Triangle), - //(c64, Laplace, hypersingular, Triangle), - //(c32, Laplace, hypersingular, Triangle), - (c64, Helmholtz, single_layer, Triangle), - (c32, Helmholtz, single_layer, Triangle), - (c64, Helmholtz, double_layer, Triangle), - (c32, Helmholtz, double_layer, Triangle), - (c64, Helmholtz, adjoint_double_layer, Triangle), - (c32, Helmholtz, adjoint_double_layer, Triangle), - (c64, Helmholtz, hypersingular, Triangle), - (c32, Helmholtz, hypersingular, Triangle), - (f64, Laplace, single_layer, Quadrilateral), - (f32, Laplace, single_layer, Quadrilateral), - //(c64, Laplace, single_layer, Quadrilateral), - //(c32, Laplace, single_layer, Quadrilateral), - (f64, Laplace, double_layer, Quadrilateral), - (f32, Laplace, double_layer, Quadrilateral), - //(c64, Laplace, double_layer, Quadrilateral), - //(c32, Laplace, double_layer, Quadrilateral), - (f64, Laplace, adjoint_double_layer, Quadrilateral), - (f32, Laplace, adjoint_double_layer, Quadrilateral), - //(c64, Laplace, adjoint_double_layer, Quadrilateral), - //(c32, Laplace, adjoint_double_layer, Quadrilateral), - (f64, Laplace, hypersingular, Quadrilateral), - (f32, Laplace, hypersingular, Quadrilateral), - //(c64, Laplace, hypersingular, Quadrilateral), - //(c32, Laplace, hypersingular, Quadrilateral), - (c64, Helmholtz, single_layer, Quadrilateral), - (c32, Helmholtz, single_layer, Quadrilateral), - (c64, Helmholtz, double_layer, Quadrilateral), - (c32, Helmholtz, double_layer, Quadrilateral), - (c64, Helmholtz, adjoint_double_layer, Quadrilateral), - (c32, Helmholtz, adjoint_double_layer, Quadrilateral), - (c64, Helmholtz, hypersingular, Quadrilateral), - (c32, Helmholtz, hypersingular, Quadrilateral) //(f64, Laplace, single_layer, Mixed), - //(f32, Laplace, single_layer, Mixed), - //(c64, Laplace, single_layer, Mixed), - //(c32, Laplace, single_layer, Mixed), - //(f64, Laplace, double_layer, Mixed), - //(f32, Laplace, double_layer, Mixed), - //(c64, Laplace, double_layer, Mixed), - //(c32, Laplace, double_layer, Mixed), - //(f64, Laplace, adjoint_double_layer, Mixed), - //(f32, Laplace, adjoint_double_layer, Mixed), - //(c64, Laplace, adjoint_double_layer, Mixed), - //(c32, Laplace, adjoint_double_layer, Mixed), - //(f64, Laplace, hypersingular, Mixed), - //(f32, Laplace, hypersingular, Mixed), - //(c64, Laplace, hypersingular, Mixed), - //(c32, Laplace, hypersingular, Mixed), - //(c64, Helmholtz, single_layer, Mixed), - //(c32, Helmholtz, single_layer, Mixed), - //(c64, Helmholtz, double_layer, Mixed), - //(c32, Helmholtz, double_layer, Mixed), - //(c64, Helmholtz, adjoint_double_layer, Mixed), - //(c32, Helmholtz, adjoint_double_layer, Mixed), - //(c64, Helmholtz, hypersingular, Mixed), - //(c32, Helmholtz, hypersingular, Mixed) - ); -} +// test_assembly!( +// (f64, Laplace, single_layer, Triangle), +// (f32, Laplace, single_layer, Triangle), +// //(c64, Laplace, single_layer, Triangle), +// //(c32, Laplace, single_layer, Triangle), +// (f64, Laplace, double_layer, Triangle), +// (f32, Laplace, double_layer, Triangle), +// //(c64, Laplace, double_layer, Triangle), +// //(c32, Laplace, double_layer, Triangle), +// (f64, Laplace, adjoint_double_layer, Triangle), +// (f32, Laplace, adjoint_double_layer, Triangle), +// //(c64, Laplace, adjoint_double_layer, Triangle), +// //(c32, Laplace, adjoint_double_layer, Triangle), +// (f64, Laplace, hypersingular, Triangle), +// (f32, Laplace, hypersingular, Triangle), +// //(c64, Laplace, hypersingular, Triangle), +// //(c32, Laplace, hypersingular, Triangle), +// (c64, Helmholtz, single_layer, Triangle), +// (c32, Helmholtz, single_layer, Triangle), +// (c64, Helmholtz, double_layer, Triangle), +// (c32, Helmholtz, double_layer, Triangle), +// (c64, Helmholtz, adjoint_double_layer, Triangle), +// (c32, Helmholtz, adjoint_double_layer, Triangle), +// (c64, Helmholtz, hypersingular, Triangle), +// (c32, Helmholtz, hypersingular, Triangle), +// (f64, Laplace, single_layer, Quadrilateral), +// (f32, Laplace, single_layer, Quadrilateral), +// //(c64, Laplace, single_layer, Quadrilateral), +// //(c32, Laplace, single_layer, Quadrilateral), +// (f64, Laplace, double_layer, Quadrilateral), +// (f32, Laplace, double_layer, Quadrilateral), +// //(c64, Laplace, double_layer, Quadrilateral), +// //(c32, Laplace, double_layer, Quadrilateral), +// (f64, Laplace, adjoint_double_layer, Quadrilateral), +// (f32, Laplace, adjoint_double_layer, Quadrilateral), +// //(c64, Laplace, adjoint_double_layer, Quadrilateral), +// //(c32, Laplace, adjoint_double_layer, Quadrilateral), +// (f64, Laplace, hypersingular, Quadrilateral), +// (f32, Laplace, hypersingular, Quadrilateral), +// //(c64, Laplace, hypersingular, Quadrilateral), +// //(c32, Laplace, hypersingular, Quadrilateral), +// (c64, Helmholtz, single_layer, Quadrilateral), +// (c32, Helmholtz, single_layer, Quadrilateral), +// (c64, Helmholtz, double_layer, Quadrilateral), +// (c32, Helmholtz, double_layer, Quadrilateral), +// (c64, Helmholtz, adjoint_double_layer, Quadrilateral), +// (c32, Helmholtz, adjoint_double_layer, Quadrilateral), +// (c64, Helmholtz, hypersingular, Quadrilateral), +// (c32, Helmholtz, hypersingular, Quadrilateral) //(f64, Laplace, single_layer, Mixed), +// //(f32, Laplace, single_layer, Mixed), +// //(c64, Laplace, single_layer, Mixed), +// //(c32, Laplace, single_layer, Mixed), +// //(f64, Laplace, double_layer, Mixed), +// //(f32, Laplace, double_layer, Mixed), +// //(c64, Laplace, double_layer, Mixed), +// //(c32, Laplace, double_layer, Mixed), +// //(f64, Laplace, adjoint_double_layer, Mixed), +// //(f32, Laplace, adjoint_double_layer, Mixed), +// //(c64, Laplace, adjoint_double_layer, Mixed), +// //(c32, Laplace, adjoint_double_layer, Mixed), +// //(f64, Laplace, hypersingular, Mixed), +// //(f32, Laplace, hypersingular, Mixed), +// //(c64, Laplace, hypersingular, Mixed), +// //(c32, Laplace, hypersingular, Mixed), +// //(c64, Helmholtz, single_layer, Mixed), +// //(c32, Helmholtz, single_layer, Mixed), +// //(c64, Helmholtz, double_layer, Mixed), +// //(c32, Helmholtz, double_layer, Mixed), +// //(c64, Helmholtz, adjoint_double_layer, Mixed), +// //(c32, Helmholtz, adjoint_double_layer, Mixed), +// //(c64, Helmholtz, hypersingular, Mixed), +// //(c32, Helmholtz, hypersingular, Mixed) +// ); +// } diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index e3ab7598..f2257d2a 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -3,14 +3,18 @@ mod assemblers; pub(crate) mod cell_pair_assemblers; pub mod integrands; -pub use assemblers::BoundaryAssembler; +pub use assemblers::{BoundaryAssembler, BoundaryAssemblerOptions}; #[cfg(test)] mod test { use super::*; + use crate::assembly::kernels::KernelEvaluator; use crate::function::SerialFunctionSpace; use crate::traits::{BoundaryAssembly, FunctionSpace}; use approx::*; + use green_kernels::laplace_3d::Laplace3dKernel; + use green_kernels::types::GreenKernelEvalType; + use integrands::SingleLayerBoundaryIntegrand; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::Continuity; use ndgrid::shapes::regular_sphere; @@ -26,7 +30,14 @@ mod test { let ndofs = space.global_size(); let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let assembler = BoundaryAssembler::::new_laplace_single_layer(); + let assembler = BoundaryAssembler::::new( + SingleLayerBoundaryIntegrand::new(), + KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value), + BoundaryAssemblerOptions::default(), + 1, + 0, + ); + // let assembler = BoundaryAssembler::::new_laplace_single_layer(); assembler.assemble_singular_into_dense(&mut matrix, &space, &space); let csr = assembler.assemble_singular_into_csr(&space, &space); @@ -52,7 +63,13 @@ mod test { let ndofs = space.global_size(); let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let assembler = BoundaryAssembler::::new_laplace_single_layer(); + let assembler = BoundaryAssembler::::new( + SingleLayerBoundaryIntegrand::new(), + KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value), + BoundaryAssemblerOptions::default(), + 1, + 0, + ); assembler.assemble_singular_into_dense(&mut matrix, &space, &space); let csr = assembler.assemble_singular_into_csr(&space, &space); @@ -81,7 +98,13 @@ mod test { let ndofs1 = space1.global_size(); let mut matrix = rlst_dynamic_array2!(f64, [ndofs1, ndofs0]); - let assembler = BoundaryAssembler::::new_laplace_single_layer(); + let assembler = BoundaryAssembler::::new( + SingleLayerBoundaryIntegrand::new(), + KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value), + BoundaryAssemblerOptions::default(), + 1, + 0, + ); assembler.assemble_singular_into_dense(&mut matrix, &space0, &space1); let csr = assembler.assemble_singular_into_csr(&space0, &space1); let indptr = csr.indptr(); diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 29148106..ac29dfc4 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -1,8 +1,4 @@ //! Boundary assemblers -pub mod adjoint_double_layer; -pub mod double_layer; -pub mod hypersingular; -pub mod single_layer; use super::cell_pair_assemblers::{ NonsingularCellPairAssembler, NonsingularCellPairAssemblerWithTestCaching, @@ -418,13 +414,14 @@ fn get_pairs_if_smallest( } /// Options for a boundary assembler +#[derive(Clone)] pub struct BoundaryAssemblerOptions { /// Number of points used in quadrature for non-singular integrals - quadrature_degrees: HashMap, + pub quadrature_degrees: HashMap, /// Quadrature degrees to be used for singular integrals - singular_quadrature_degrees: HashMap<(ReferenceCellType, ReferenceCellType), usize>, + pub singular_quadrature_degrees: HashMap<(ReferenceCellType, ReferenceCellType), usize>, /// Maximum size of each batch of cells to send to an assembly function - batch_size: usize, + pub batch_size: usize, } impl Default for BoundaryAssemblerOptions { @@ -443,6 +440,31 @@ impl Default for BoundaryAssemblerOptions { } } +impl BoundaryAssemblerOptions { + /// Set the regular quadrature order. + pub fn set_regular_quadrature_degree(&mut self, cell_type: ReferenceCellType, npoints: usize) { + self.quadrature_degrees + .entry(cell_type) + .and_modify(|x| *x = npoints); + } + + /// Set the singular quadrature order. + pub fn set_singular_quadrature_degree( + &mut self, + cell_type: (ReferenceCellType, ReferenceCellType), + npoints: usize, + ) { + self.singular_quadrature_degrees + .entry(cell_type) + .and_modify(|x| *x = npoints); + } + + /// Set the batch size. + pub fn set_batch_size(&mut self, batch_size: usize) { + self.batch_size = batch_size; + } +} + /// Boundary assembler /// /// Assembles operators by processing batches of cells in parallel @@ -473,11 +495,17 @@ impl< > BoundaryAssembler { /// Create new - fn new(integrand: Integrand, kernel: Kernel, deriv_size: usize, table_derivs: usize) -> Self { + pub fn new( + integrand: Integrand, + kernel: Kernel, + options: BoundaryAssemblerOptions, + deriv_size: usize, + table_derivs: usize, + ) -> Self { Self { integrand, kernel, - options: BoundaryAssemblerOptions::default(), + options, deriv_size, table_derivs, } diff --git a/src/assembly/boundary/assemblers/adjoint_double_layer.rs b/src/assembly/boundary/assemblers/adjoint_double_layer.rs deleted file mode 100644 index f1e01018..00000000 --- a/src/assembly/boundary/assemblers/adjoint_double_layer.rs +++ /dev/null @@ -1,46 +0,0 @@ -//! Adjoint double layer assemblers -use super::BoundaryAssembler; -use crate::assembly::{ - boundary::integrands::AdjointDoubleLayerBoundaryIntegrand, common::GreenKernelEvalType, - kernels::KernelEvaluator, -}; -use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar}; - -impl> - BoundaryAssembler, KernelEvaluator> -{ - /// Create a new adjoint double layer assembler - pub fn new_adjoint_double_layer(kernel: KernelEvaluator) -> Self { - Self::new(AdjointDoubleLayerBoundaryIntegrand::new(), kernel, 4, 0) - } -} -impl - BoundaryAssembler< - T, - AdjointDoubleLayerBoundaryIntegrand, - KernelEvaluator>, - > -{ - /// Create a new Laplace adjoint double layer assembler - pub fn new_laplace_adjoint_double_layer() -> Self { - Self::new_adjoint_double_layer(KernelEvaluator::new_laplace( - GreenKernelEvalType::ValueDeriv, - )) - } -} -impl + MatrixInverse> - BoundaryAssembler< - T, - AdjointDoubleLayerBoundaryIntegrand, - KernelEvaluator>, - > -{ - /// Create a new Helmholtz adjoint double layer assembler - pub fn new_helmholtz_adjoint_double_layer(wavenumber: T::Real) -> Self { - Self::new_adjoint_double_layer(KernelEvaluator::new_helmholtz( - wavenumber, - GreenKernelEvalType::ValueDeriv, - )) - } -} diff --git a/src/assembly/boundary/assemblers/double_layer.rs b/src/assembly/boundary/assemblers/double_layer.rs deleted file mode 100644 index f837ea5e..00000000 --- a/src/assembly/boundary/assemblers/double_layer.rs +++ /dev/null @@ -1,38 +0,0 @@ -//! Double layer assemblers -use super::BoundaryAssembler; -use crate::assembly::{ - boundary::integrands::DoubleLayerBoundaryIntegrand, common::GreenKernelEvalType, - kernels::KernelEvaluator, -}; -use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar}; - -impl> - BoundaryAssembler, KernelEvaluator> -{ - /// Create a new double layer assembler - pub fn new_double_layer(kernel: KernelEvaluator) -> Self { - Self::new(DoubleLayerBoundaryIntegrand::new(), kernel, 4, 0) - } -} -impl - BoundaryAssembler, KernelEvaluator>> -{ - /// Create a new Laplace double layer assembler - pub fn new_laplace_double_layer() -> Self { - Self::new_double_layer(KernelEvaluator::new_laplace( - GreenKernelEvalType::ValueDeriv, - )) - } -} -impl + MatrixInverse> - BoundaryAssembler, KernelEvaluator>> -{ - /// Create a new Helmholtz double layer assembler - pub fn new_helmholtz_double_layer(wavenumber: T::Real) -> Self { - Self::new_double_layer(KernelEvaluator::new_helmholtz( - wavenumber, - GreenKernelEvalType::ValueDeriv, - )) - } -} diff --git a/src/assembly/boundary/assemblers/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs deleted file mode 100644 index d53684de..00000000 --- a/src/assembly/boundary/assemblers/hypersingular.rs +++ /dev/null @@ -1,60 +0,0 @@ -//! Hypersingular assemblers -use super::BoundaryAssembler; -use crate::assembly::{ - boundary::integrands::{ - BoundaryIntegrandScalarProduct, BoundaryIntegrandSum, - HypersingularCurlCurlBoundaryIntegrand, HypersingularNormalNormalBoundaryIntegrand, - }, - common::GreenKernelEvalType, - kernels::KernelEvaluator, -}; -use crate::traits::BoundaryIntegrand; -use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar}; - -type HelmholtzIntegrand = BoundaryIntegrandSum< - T, - HypersingularCurlCurlBoundaryIntegrand, - BoundaryIntegrandScalarProduct>, ->; - -impl, I: BoundaryIntegrand> - BoundaryAssembler> -{ - /// Create a new adjoint double layer assembler - pub fn new_hypersingular(integrand: I, kernel: KernelEvaluator) -> Self { - Self::new(integrand, kernel, 4, 1) - } -} -impl - BoundaryAssembler< - T, - HypersingularCurlCurlBoundaryIntegrand, - KernelEvaluator>, - > -{ - /// Create a new Laplace adjoint double layer assembler - pub fn new_laplace_hypersingular() -> Self { - Self::new_hypersingular( - HypersingularCurlCurlBoundaryIntegrand::new(), - KernelEvaluator::new_laplace(GreenKernelEvalType::ValueDeriv), - ) - } -} -impl + MatrixInverse> - BoundaryAssembler, KernelEvaluator>> -{ - /// Create a new Helmholtz adjoint double layer assembler - pub fn new_helmholtz_hypersingular(wavenumber: T::Real) -> Self { - Self::new_hypersingular( - BoundaryIntegrandSum::new( - HypersingularCurlCurlBoundaryIntegrand::new(), - BoundaryIntegrandScalarProduct::new( - num::cast::(-wavenumber.powi(2)).unwrap(), - HypersingularNormalNormalBoundaryIntegrand::new(), - ), - ), - KernelEvaluator::new_helmholtz(wavenumber, GreenKernelEvalType::ValueDeriv), - ) - } -} diff --git a/src/assembly/boundary/assemblers/single_layer.rs b/src/assembly/boundary/assemblers/single_layer.rs deleted file mode 100644 index 90045f6d..00000000 --- a/src/assembly/boundary/assemblers/single_layer.rs +++ /dev/null @@ -1,36 +0,0 @@ -//! Single layer assemblers -use super::BoundaryAssembler; -use crate::assembly::{ - boundary::integrands::SingleLayerBoundaryIntegrand, common::GreenKernelEvalType, - kernels::KernelEvaluator, -}; -use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar}; - -impl> - BoundaryAssembler, KernelEvaluator> -{ - /// Create a new single layer assembler - pub fn new_single_layer(kernel: KernelEvaluator) -> Self { - Self::new(SingleLayerBoundaryIntegrand::new(), kernel, 1, 0) - } -} -impl - BoundaryAssembler, KernelEvaluator>> -{ - /// Create a new Laplace single layer assembler - pub fn new_laplace_single_layer() -> Self { - Self::new_single_layer(KernelEvaluator::new_laplace(GreenKernelEvalType::Value)) - } -} -impl + MatrixInverse> - BoundaryAssembler, KernelEvaluator>> -{ - /// Create a new Helmholtz single layer assembler - pub fn new_helmholtz_single_layer(wavenumber: T::Real) -> Self { - Self::new_single_layer(KernelEvaluator::new_helmholtz( - wavenumber, - GreenKernelEvalType::Value, - )) - } -} diff --git a/src/bindings.rs b/src/bindings.rs index ebac62ef..1737fdbe 100644 --- a/src/bindings.rs +++ b/src/bindings.rs @@ -18,7 +18,7 @@ impl DType { 0 => Some(DType::F32), 1 => Some(DType::F64), 2 => Some(DType::C32), - 3 => Some(DType::C64), + _ => None, } } @@ -4334,1094 +4334,1094 @@ pub mod boundary_assembly { ) -> u8 { (*assembler).dtype as u8 } +} +// unsafe fn laplace_boundary_assembler_new_internal( +// operator: BoundaryOperator, +// dtype: DType, +// ) -> *const BoundaryAssemblerWrapper { +// Box::into_raw(Box::new(BoundaryAssemblerWrapper { +// assembler: match operator { +// BoundaryOperator::SingleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_laplace_single_layer(), +// )) as *const c_void, +// BoundaryOperator::DoubleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_laplace_double_layer(), +// )) as *const c_void, +// BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_laplace_adjoint_double_layer(), +// )) as *const c_void, +// BoundaryOperator::Hypersingular => Box::into_raw(Box::new( +// BoundaryAssembler::::new_laplace_hypersingular(), +// )) as *const c_void, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// itype: operator, +// ktype: KernelType::Laplace, +// dtype, +// })) +// } - unsafe fn laplace_boundary_assembler_new_internal( - operator: BoundaryOperator, - dtype: DType, - ) -> *const BoundaryAssemblerWrapper { - Box::into_raw(Box::new(BoundaryAssemblerWrapper { - assembler: match operator { - BoundaryOperator::SingleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_laplace_single_layer(), - )) as *const c_void, - BoundaryOperator::DoubleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_laplace_double_layer(), - )) as *const c_void, - BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_laplace_adjoint_double_layer(), - )) as *const c_void, - BoundaryOperator::Hypersingular => Box::into_raw(Box::new( - BoundaryAssembler::::new_laplace_hypersingular(), - )) as *const c_void, - _ => { - panic!("Invalid operator"); - } - }, - itype: operator, - ktype: KernelType::Laplace, - dtype, - })) - } +// #[no_mangle] +// pub unsafe extern "C" fn laplace_boundary_assembler_new( +// operator: u8, +// dtype: u8, +// ) -> *const BoundaryAssemblerWrapper { +// let operator = BoundaryOperator::from(operator).unwrap(); +// let dtype = DType::from(dtype).unwrap(); +// match dtype { +// DType::F32 => laplace_boundary_assembler_new_internal::(operator, dtype), +// DType::F64 => laplace_boundary_assembler_new_internal::(operator, dtype), +// _ => { +// panic!("Invalid data type"); +// } +// } +// } - #[no_mangle] - pub unsafe extern "C" fn laplace_boundary_assembler_new( - operator: u8, - dtype: u8, - ) -> *const BoundaryAssemblerWrapper { - let operator = BoundaryOperator::from(operator).unwrap(); - let dtype = DType::from(dtype).unwrap(); - match dtype { - DType::F32 => laplace_boundary_assembler_new_internal::(operator, dtype), - DType::F64 => laplace_boundary_assembler_new_internal::(operator, dtype), - _ => { - panic!("Invalid data type"); - } - } - } +// unsafe fn helmholtz_boundary_assembler_new_internal< +// T: RlstScalar + MatrixInverse, +// >( +// wavenumber: T::Real, +// operator: BoundaryOperator, +// dtype: DType, +// ) -> *const BoundaryAssemblerWrapper { +// Box::into_raw(Box::new(BoundaryAssemblerWrapper { +// assembler: match operator { +// BoundaryOperator::SingleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_helmholtz_single_layer(wavenumber), +// )) as *const c_void, +// BoundaryOperator::DoubleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_helmholtz_double_layer(wavenumber), +// )) as *const c_void, +// BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new( +// BoundaryAssembler::::new_helmholtz_adjoint_double_layer(wavenumber), +// )) as *const c_void, +// BoundaryOperator::Hypersingular => Box::into_raw(Box::new( +// BoundaryAssembler::::new_helmholtz_hypersingular(wavenumber), +// )) as *const c_void, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// itype: operator, +// ktype: KernelType::Helmholtz, +// dtype, +// })) +// } - unsafe fn helmholtz_boundary_assembler_new_internal< - T: RlstScalar + MatrixInverse, - >( - wavenumber: T::Real, - operator: BoundaryOperator, - dtype: DType, - ) -> *const BoundaryAssemblerWrapper { - Box::into_raw(Box::new(BoundaryAssemblerWrapper { - assembler: match operator { - BoundaryOperator::SingleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_helmholtz_single_layer(wavenumber), - )) as *const c_void, - BoundaryOperator::DoubleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_helmholtz_double_layer(wavenumber), - )) as *const c_void, - BoundaryOperator::AdjointDoubleLayer => Box::into_raw(Box::new( - BoundaryAssembler::::new_helmholtz_adjoint_double_layer(wavenumber), - )) as *const c_void, - BoundaryOperator::Hypersingular => Box::into_raw(Box::new( - BoundaryAssembler::::new_helmholtz_hypersingular(wavenumber), - )) as *const c_void, - _ => { - panic!("Invalid operator"); - } - }, - itype: operator, - ktype: KernelType::Helmholtz, - dtype, - })) - } +// #[no_mangle] +// pub unsafe extern "C" fn helmholtz_boundary_assembler_new( +// wavenumber: *const c_void, +// operator: u8, +// dtype: u8, +// ) -> *const BoundaryAssemblerWrapper { +// let operator = BoundaryOperator::from(operator).unwrap(); +// let dtype = DType::from(dtype).unwrap(); +// match dtype { +// DType::C32 => helmholtz_boundary_assembler_new_internal::( +// *(wavenumber as *const f32), +// operator, +// dtype, +// ), +// DType::C64 => helmholtz_boundary_assembler_new_internal::( +// *(wavenumber as *const f64), +// operator, +// dtype, +// ), +// _ => { +// panic!("Invalid data type"); +// } +// } +// } +// } - #[no_mangle] - pub unsafe extern "C" fn helmholtz_boundary_assembler_new( - wavenumber: *const c_void, - operator: u8, - dtype: u8, - ) -> *const BoundaryAssemblerWrapper { - let operator = BoundaryOperator::from(operator).unwrap(); - let dtype = DType::from(dtype).unwrap(); - match dtype { - DType::C32 => helmholtz_boundary_assembler_new_internal::( - *(wavenumber as *const f32), - operator, - dtype, - ), - DType::C64 => helmholtz_boundary_assembler_new_internal::( - *(wavenumber as *const f64), - operator, - dtype, - ), - _ => { - panic!("Invalid data type"); - } - } - } -} +// pub mod potential_assembly { +// use super::boundary_assembly::KernelType; +// use super::function::{extract_space, FunctionSpaceWrapper, GridType, SpaceType}; +// use super::DType; +// use crate::{ +// assembly::kernels::KernelEvaluator, +// assembly::potential::integrands::{ +// DoubleLayerPotentialIntegrand, SingleLayerPotentialIntegrand, +// }, +// assembly::potential::PotentialAssembler, +// function::SerialFunctionSpace, +// traits::{ +// FunctionSpace, KernelEvaluator as KernelEvaluatorTrait, PotentialAssembly, +// PotentialIntegrand, +// }, +// }; +// use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel}; +// use ndelement::{ciarlet::CiarletElement, types::ReferenceCellType}; +// use ndgrid::SingleElementGrid; +// use rlst::{ +// c32, c64, rlst_array_from_slice2, rlst_array_from_slice_mut2, MatrixInverse, RlstScalar, +// }; +// use std::ffi::c_void; +// use std::slice::{from_raw_parts, from_raw_parts_mut}; -pub mod potential_assembly { - use super::boundary_assembly::KernelType; - use super::function::{extract_space, FunctionSpaceWrapper, GridType, SpaceType}; - use super::DType; - use crate::{ - assembly::kernels::KernelEvaluator, - assembly::potential::integrands::{ - DoubleLayerPotentialIntegrand, SingleLayerPotentialIntegrand, - }, - assembly::potential::PotentialAssembler, - function::SerialFunctionSpace, - traits::{ - FunctionSpace, KernelEvaluator as KernelEvaluatorTrait, PotentialAssembly, - PotentialIntegrand, - }, - }; - use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel}; - use ndelement::{ciarlet::CiarletElement, types::ReferenceCellType}; - use ndgrid::SingleElementGrid; - use rlst::{ - c32, c64, rlst_array_from_slice2, rlst_array_from_slice_mut2, MatrixInverse, RlstScalar, - }; - use std::ffi::c_void; - use std::slice::{from_raw_parts, from_raw_parts_mut}; +// #[derive(Debug, PartialEq, Clone, Copy)] +// #[repr(u8)] +// pub enum PotentialOperator { +// SingleLayer = 0, +// DoubleLayer = 1, +// ElectricField = 2, +// MagneticField = 3, +// } - #[derive(Debug, PartialEq, Clone, Copy)] - #[repr(u8)] - pub enum PotentialOperator { - SingleLayer = 0, - DoubleLayer = 1, - ElectricField = 2, - MagneticField = 3, - } +// impl PotentialOperator { +// fn from(value: u8) -> Option { +// match value { +// 0 => Some(PotentialOperator::SingleLayer), +// 1 => Some(PotentialOperator::DoubleLayer), +// 2 => Some(PotentialOperator::ElectricField), +// 3 => Some(PotentialOperator::MagneticField), +// _ => None, +// } +// } +// } - impl PotentialOperator { - fn from(value: u8) -> Option { - match value { - 0 => Some(PotentialOperator::SingleLayer), - 1 => Some(PotentialOperator::DoubleLayer), - 2 => Some(PotentialOperator::ElectricField), - 3 => Some(PotentialOperator::MagneticField), - _ => None, - } - } - } +// #[repr(C)] +// pub struct PotentialAssemblerWrapper { +// pub assembler: *const c_void, +// pub itype: PotentialOperator, +// pub ktype: KernelType, +// pub dtype: DType, +// } +// impl Drop for PotentialAssemblerWrapper { +// fn drop(&mut self) { +// let Self { +// assembler, +// itype, +// ktype, +// dtype, +// } = self; +// match ktype { +// KernelType::Laplace => match itype { +// PotentialOperator::SingleLayer => match dtype { +// DType::F32 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// DType::F64 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match dtype { +// DType::F32 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// DType::F64 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match itype { +// PotentialOperator::SingleLayer => match dtype { +// DType::C32 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// DType::C64 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match dtype { +// DType::C32 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// DType::C64 => drop(unsafe { +// Box::from_raw( +// *assembler +// as *mut PotentialAssembler< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >, +// ) +// }), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// } +// } - #[repr(C)] - pub struct PotentialAssemblerWrapper { - pub assembler: *const c_void, - pub itype: PotentialOperator, - pub ktype: KernelType, - pub dtype: DType, - } - impl Drop for PotentialAssemblerWrapper { - fn drop(&mut self) { - let Self { - assembler, - itype, - ktype, - dtype, - } = self; - match ktype { - KernelType::Laplace => match itype { - PotentialOperator::SingleLayer => match dtype { - DType::F32 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - DType::F64 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match dtype { - DType::F32 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - DType::F64 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match itype { - PotentialOperator::SingleLayer => match dtype { - DType::C32 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - DType::C64 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match dtype { - DType::C32 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - DType::C64 => drop(unsafe { - Box::from_raw( - *assembler - as *mut PotentialAssembler< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >, - ) - }), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - } - } +// #[no_mangle] +// pub unsafe extern "C" fn free_potential_assembler(a: *mut PotentialAssemblerWrapper) { +// assert!(!a.is_null()); +// unsafe { drop(Box::from_raw(a)) } +// } - #[no_mangle] - pub unsafe extern "C" fn free_potential_assembler(a: *mut PotentialAssemblerWrapper) { - assert!(!a.is_null()); - unsafe { drop(Box::from_raw(a)) } - } +// pub(crate) unsafe fn extract_potential_assembler< +// T: RlstScalar + MatrixInverse, +// Integrand: PotentialIntegrand, +// Kernel: KernelEvaluatorTrait, +// >( +// assembler: *const PotentialAssemblerWrapper, +// ) -> *const PotentialAssembler { +// (*assembler).assembler as *const PotentialAssembler +// } - pub(crate) unsafe fn extract_potential_assembler< - T: RlstScalar + MatrixInverse, - Integrand: PotentialIntegrand, - Kernel: KernelEvaluatorTrait, - >( - assembler: *const PotentialAssemblerWrapper, - ) -> *const PotentialAssembler { - (*assembler).assembler as *const PotentialAssembler - } +// pub(crate) unsafe fn extract_potential_assembler_mut< +// T: RlstScalar + MatrixInverse, +// Integrand: PotentialIntegrand, +// Kernel: KernelEvaluatorTrait, +// >( +// assembler: *const PotentialAssemblerWrapper, +// ) -> *mut PotentialAssembler { +// (*assembler).assembler as *mut PotentialAssembler +// } - pub(crate) unsafe fn extract_potential_assembler_mut< - T: RlstScalar + MatrixInverse, - Integrand: PotentialIntegrand, - Kernel: KernelEvaluatorTrait, - >( - assembler: *const PotentialAssemblerWrapper, - ) -> *mut PotentialAssembler { - (*assembler).assembler as *mut PotentialAssembler - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_has_quadrature_degree( +// assembler: *mut PotentialAssemblerWrapper, +// cell: u8, +// ) -> bool { +// let cell = ReferenceCellType::from(cell).unwrap(); +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// .is_some() +// } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_has_quadrature_degree( - assembler: *mut PotentialAssemblerWrapper, - cell: u8, - ) -> bool { - let cell = ReferenceCellType::from(cell).unwrap(); - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::F64 => (*extract_potential_assembler::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::F64 => (*extract_potential_assembler::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::C64 => (*extract_potential_assembler::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::C64 => (*extract_potential_assembler::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - .is_some() - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_set_quadrature_degree( +// assembler: *mut PotentialAssemblerWrapper, +// cell: u8, +// degree: usize, +// ) { +// let cell = ReferenceCellType::from(cell).unwrap(); +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler_mut::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// DType::F64 => (*extract_potential_assembler_mut::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler_mut::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// DType::F64 => (*extract_potential_assembler_mut::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler_mut::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// DType::C64 => (*extract_potential_assembler_mut::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler_mut::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// DType::C64 => (*extract_potential_assembler_mut::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_quadrature_degree(cell, degree), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_set_quadrature_degree( - assembler: *mut PotentialAssemblerWrapper, - cell: u8, - degree: usize, - ) { - let cell = ReferenceCellType::from(cell).unwrap(); - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler_mut::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - DType::F64 => (*extract_potential_assembler_mut::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler_mut::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - DType::F64 => (*extract_potential_assembler_mut::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler_mut::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - DType::C64 => (*extract_potential_assembler_mut::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler_mut::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - DType::C64 => (*extract_potential_assembler_mut::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_quadrature_degree(cell, degree), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - } - - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_quadrature_degree( - assembler: *mut PotentialAssemblerWrapper, - cell: u8, - ) -> usize { - let cell = ReferenceCellType::from(cell).unwrap(); - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::F64 => (*extract_potential_assembler::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::F64 => (*extract_potential_assembler::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::C64 => (*extract_potential_assembler::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - DType::C64 => (*extract_potential_assembler::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .quadrature_degree(cell), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - .unwrap() - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_quadrature_degree( +// assembler: *mut PotentialAssemblerWrapper, +// cell: u8, +// ) -> usize { +// let cell = ReferenceCellType::from(cell).unwrap(); +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .quadrature_degree(cell), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// .unwrap() +// } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_set_batch_size( - assembler: *mut PotentialAssemblerWrapper, - batch_size: usize, - ) { - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler_mut::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - DType::F64 => (*extract_potential_assembler_mut::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler_mut::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - DType::F64 => (*extract_potential_assembler_mut::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler_mut::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - DType::C64 => (*extract_potential_assembler_mut::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler_mut::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - DType::C64 => (*extract_potential_assembler_mut::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .set_batch_size(batch_size), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_set_batch_size( +// assembler: *mut PotentialAssemblerWrapper, +// batch_size: usize, +// ) { +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler_mut::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// DType::F64 => (*extract_potential_assembler_mut::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler_mut::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// DType::F64 => (*extract_potential_assembler_mut::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler_mut::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// DType::C64 => (*extract_potential_assembler_mut::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler_mut::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// DType::C64 => (*extract_potential_assembler_mut::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .set_batch_size(batch_size), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_batch_size( - assembler: *mut PotentialAssemblerWrapper, - ) -> usize { - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - DType::F64 => (*extract_potential_assembler::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => (*extract_potential_assembler::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - DType::F64 => (*extract_potential_assembler::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - DType::C64 => (*extract_potential_assembler::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => (*extract_potential_assembler::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - DType::C64 => (*extract_potential_assembler::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - >(assembler)) - .batch_size(), - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_batch_size( +// assembler: *mut PotentialAssemblerWrapper, +// ) -> usize { +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => (*extract_potential_assembler::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// DType::F64 => (*extract_potential_assembler::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => (*extract_potential_assembler::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// DType::C64 => (*extract_potential_assembler::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// >(assembler)) +// .batch_size(), +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// } - unsafe fn potential_assembler_assemble_into_dense_internal_real< - T: RlstScalar + MatrixInverse, - Integrand: PotentialIntegrand, - Kernel: KernelEvaluatorTrait, - Space: FunctionSpace + Sync, - >( - assembler: *mut PotentialAssemblerWrapper, - output: *mut c_void, - space: *const FunctionSpaceWrapper, - points: *const c_void, - npts: usize, - ) { - let points = rlst_array_from_slice2!( - from_raw_parts(points as *const T::Real, 3 * npts), - [3, npts] - ); - let dim = (*extract_space::(space)).global_size(); - let mut output = rlst_array_from_slice_mut2!( - from_raw_parts_mut(output as *mut T, npts * dim), - [npts, dim] - ); +// unsafe fn potential_assembler_assemble_into_dense_internal_real< +// T: RlstScalar + MatrixInverse, +// Integrand: PotentialIntegrand, +// Kernel: KernelEvaluatorTrait, +// Space: FunctionSpace + Sync, +// >( +// assembler: *mut PotentialAssemblerWrapper, +// output: *mut c_void, +// space: *const FunctionSpaceWrapper, +// points: *const c_void, +// npts: usize, +// ) { +// let points = rlst_array_from_slice2!( +// from_raw_parts(points as *const T::Real, 3 * npts), +// [3, npts] +// ); +// let dim = (*extract_space::(space)).global_size(); +// let mut output = rlst_array_from_slice_mut2!( +// from_raw_parts_mut(output as *mut T, npts * dim), +// [npts, dim] +// ); - (*extract_potential_assembler::(assembler)).assemble_into_dense( - &mut output, - &*extract_space::(space), - &points, - ) - } - unsafe fn potential_assembler_assemble_into_dense_internal_complex< - T: RlstScalar + MatrixInverse, - Integrand: PotentialIntegrand, - Kernel: KernelEvaluatorTrait, - Space: FunctionSpace + Sync, - >( - assembler: *mut PotentialAssemblerWrapper, - output: *mut c_void, - space: *const FunctionSpaceWrapper, - points: *const c_void, - npts: usize, - ) { - let points = rlst_array_from_slice2!( - from_raw_parts(points as *const T::Real, 3 * npts), - [3, npts] - ); - let dim = (*extract_space::(space)).global_size(); - let mut output = rlst_array_from_slice_mut2!( - from_raw_parts_mut(output as *mut T, npts * dim), - [npts, dim] - ); +// (*extract_potential_assembler::(assembler)).assemble_into_dense( +// &mut output, +// &*extract_space::(space), +// &points, +// ) +// } +// unsafe fn potential_assembler_assemble_into_dense_internal_complex< +// T: RlstScalar + MatrixInverse, +// Integrand: PotentialIntegrand, +// Kernel: KernelEvaluatorTrait, +// Space: FunctionSpace + Sync, +// >( +// assembler: *mut PotentialAssemblerWrapper, +// output: *mut c_void, +// space: *const FunctionSpaceWrapper, +// points: *const c_void, +// npts: usize, +// ) { +// let points = rlst_array_from_slice2!( +// from_raw_parts(points as *const T::Real, 3 * npts), +// [3, npts] +// ); +// let dim = (*extract_space::(space)).global_size(); +// let mut output = rlst_array_from_slice_mut2!( +// from_raw_parts_mut(output as *mut T, npts * dim), +// [npts, dim] +// ); - (*extract_potential_assembler::(assembler)).assemble_into_dense( - &mut output, - &*extract_space::(space), - &points, - ) - } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_assemble_into_dense( - assembler: *mut PotentialAssemblerWrapper, - output: *mut c_void, - space: *const FunctionSpaceWrapper, - points: *const c_void, - npts: usize, - ) { - match (*assembler).dtype { - DType::F32 => { - assert_eq!((*space).dtype, DType::F32); - } - DType::F64 => { - assert_eq!((*space).dtype, DType::F64); - } - DType::C32 => { - assert_eq!((*space).dtype, DType::F32); - } - DType::C64 => { - assert_eq!((*space).dtype, DType::F64); - } - } - match (*assembler).ktype { - KernelType::Laplace => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::F32 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_real::< - f32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - f32, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - DType::F64 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_real::< - f64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - f64, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::F32 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_real::< - f32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - f32, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - DType::F64 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_real::< - f64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - f64, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - KernelType::Helmholtz => match (*assembler).itype { - PotentialOperator::SingleLayer => match (*assembler).dtype { - DType::C32 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_complex::< - c32, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - c32, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - DType::C64 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_complex::< - c64, - SingleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - c64, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - _ => { - panic!("Invalid data type"); - } - }, - PotentialOperator::DoubleLayer => match (*assembler).dtype { - DType::C32 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_complex::< - c32, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - c32, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - DType::C64 => match (*space).stype { - SpaceType::SerialFunctionSpace => match (*space).gtype { - GridType::SerialSingleElementGrid => { - potential_assembler_assemble_into_dense_internal_complex::< - c64, - DoubleLayerPotentialIntegrand, - KernelEvaluator>, - SerialFunctionSpace< - c64, - SingleElementGrid>, - >, - >( - assembler, output, space, points, npts - ); - } - }, - }, - _ => { - panic!("Invalid data type"); - } - }, - _ => { - panic!("Invalid operator"); - } - }, - } - } +// (*extract_potential_assembler::(assembler)).assemble_into_dense( +// &mut output, +// &*extract_space::(space), +// &points, +// ) +// } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_assemble_into_dense( +// assembler: *mut PotentialAssemblerWrapper, +// output: *mut c_void, +// space: *const FunctionSpaceWrapper, +// points: *const c_void, +// npts: usize, +// ) { +// match (*assembler).dtype { +// DType::F32 => { +// assert_eq!((*space).dtype, DType::F32); +// } +// DType::F64 => { +// assert_eq!((*space).dtype, DType::F64); +// } +// DType::C32 => { +// assert_eq!((*space).dtype, DType::F32); +// } +// DType::C64 => { +// assert_eq!((*space).dtype, DType::F64); +// } +// } +// match (*assembler).ktype { +// KernelType::Laplace => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::F32 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_real::< +// f32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// f32, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// DType::F64 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_real::< +// f64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// f64, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::F32 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_real::< +// f32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// f32, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// DType::F64 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_real::< +// f64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// f64, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// KernelType::Helmholtz => match (*assembler).itype { +// PotentialOperator::SingleLayer => match (*assembler).dtype { +// DType::C32 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_complex::< +// c32, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// c32, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// DType::C64 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_complex::< +// c64, +// SingleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// c64, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// _ => { +// panic!("Invalid data type"); +// } +// }, +// PotentialOperator::DoubleLayer => match (*assembler).dtype { +// DType::C32 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_complex::< +// c32, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// c32, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// DType::C64 => match (*space).stype { +// SpaceType::SerialFunctionSpace => match (*space).gtype { +// GridType::SerialSingleElementGrid => { +// potential_assembler_assemble_into_dense_internal_complex::< +// c64, +// DoubleLayerPotentialIntegrand, +// KernelEvaluator>, +// SerialFunctionSpace< +// c64, +// SingleElementGrid>, +// >, +// >( +// assembler, output, space, points, npts +// ); +// } +// }, +// }, +// _ => { +// panic!("Invalid data type"); +// } +// }, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// } +// } - #[no_mangle] - pub unsafe extern "C" fn potential_assembler_dtype( - assembler: *mut PotentialAssemblerWrapper, - ) -> u8 { - (*assembler).dtype as u8 - } +// #[no_mangle] +// pub unsafe extern "C" fn potential_assembler_dtype( +// assembler: *mut PotentialAssemblerWrapper, +// ) -> u8 { +// (*assembler).dtype as u8 +// } - unsafe fn laplace_potential_assembler_new_internal( - operator: PotentialOperator, - dtype: DType, - ) -> *const PotentialAssemblerWrapper { - Box::into_raw(Box::new(PotentialAssemblerWrapper { - assembler: match operator { - PotentialOperator::SingleLayer => { - Box::into_raw(Box::new( - PotentialAssembler::::new_laplace_single_layer(), - )) as *const c_void - } - PotentialOperator::DoubleLayer => { - Box::into_raw(Box::new( - PotentialAssembler::::new_laplace_double_layer(), - )) as *const c_void - } - _ => { - panic!("Invalid operator"); - } - }, - itype: operator, - ktype: KernelType::Laplace, - dtype, - })) - } +// unsafe fn laplace_potential_assembler_new_internal( +// operator: PotentialOperator, +// dtype: DType, +// ) -> *const PotentialAssemblerWrapper { +// Box::into_raw(Box::new(PotentialAssemblerWrapper { +// assembler: match operator { +// PotentialOperator::SingleLayer => { +// Box::into_raw(Box::new( +// PotentialAssembler::::new_laplace_single_layer(), +// )) as *const c_void +// } +// PotentialOperator::DoubleLayer => { +// Box::into_raw(Box::new( +// PotentialAssembler::::new_laplace_double_layer(), +// )) as *const c_void +// } +// _ => { +// panic!("Invalid operator"); +// } +// }, +// itype: operator, +// ktype: KernelType::Laplace, +// dtype, +// })) +// } - #[no_mangle] - pub unsafe extern "C" fn laplace_potential_assembler_new( - operator: u8, - dtype: u8, - ) -> *const PotentialAssemblerWrapper { - let operator = PotentialOperator::from(operator).unwrap(); - let dtype = DType::from(dtype).unwrap(); - match dtype { - DType::F32 => laplace_potential_assembler_new_internal::(operator, dtype), - DType::F64 => laplace_potential_assembler_new_internal::(operator, dtype), - _ => { - panic!("Invalid data type"); - } - } - } +// #[no_mangle] +// pub unsafe extern "C" fn laplace_potential_assembler_new( +// operator: u8, +// dtype: u8, +// ) -> *const PotentialAssemblerWrapper { +// let operator = PotentialOperator::from(operator).unwrap(); +// let dtype = DType::from(dtype).unwrap(); +// match dtype { +// DType::F32 => laplace_potential_assembler_new_internal::(operator, dtype), +// DType::F64 => laplace_potential_assembler_new_internal::(operator, dtype), +// _ => { +// panic!("Invalid data type"); +// } +// } +// } - unsafe fn helmholtz_potential_assembler_new_internal< - T: RlstScalar + MatrixInverse, - >( - wavenumber: T::Real, - operator: PotentialOperator, - dtype: DType, - ) -> *const PotentialAssemblerWrapper { - Box::into_raw(Box::new(PotentialAssemblerWrapper { - assembler: match operator { - PotentialOperator::SingleLayer => Box::into_raw(Box::new( - PotentialAssembler::::new_helmholtz_single_layer(wavenumber), - )) as *const c_void, - PotentialOperator::DoubleLayer => Box::into_raw(Box::new( - PotentialAssembler::::new_helmholtz_double_layer(wavenumber), - )) as *const c_void, - _ => { - panic!("Invalid operator"); - } - }, - itype: operator, - ktype: KernelType::Helmholtz, - dtype, - })) - } +// unsafe fn helmholtz_potential_assembler_new_internal< +// T: RlstScalar + MatrixInverse, +// >( +// wavenumber: T::Real, +// operator: PotentialOperator, +// dtype: DType, +// ) -> *const PotentialAssemblerWrapper { +// Box::into_raw(Box::new(PotentialAssemblerWrapper { +// assembler: match operator { +// PotentialOperator::SingleLayer => Box::into_raw(Box::new( +// PotentialAssembler::::new_helmholtz_single_layer(wavenumber), +// )) as *const c_void, +// PotentialOperator::DoubleLayer => Box::into_raw(Box::new( +// PotentialAssembler::::new_helmholtz_double_layer(wavenumber), +// )) as *const c_void, +// _ => { +// panic!("Invalid operator"); +// } +// }, +// itype: operator, +// ktype: KernelType::Helmholtz, +// dtype, +// })) +// } - #[no_mangle] - pub unsafe extern "C" fn helmholtz_potential_assembler_new( - wavenumber: *const c_void, - operator: u8, - dtype: u8, - ) -> *const PotentialAssemblerWrapper { - let operator = PotentialOperator::from(operator).unwrap(); - let dtype = DType::from(dtype).unwrap(); - match dtype { - DType::C32 => helmholtz_potential_assembler_new_internal::( - *(wavenumber as *const f32), - operator, - dtype, - ), - DType::C64 => helmholtz_potential_assembler_new_internal::( - *(wavenumber as *const f64), - operator, - dtype, - ), - _ => { - panic!("Invalid data type"); - } - } - } -} +// #[no_mangle] +// pub unsafe extern "C" fn helmholtz_potential_assembler_new( +// wavenumber: *const c_void, +// operator: u8, +// dtype: u8, +// ) -> *const PotentialAssemblerWrapper { +// let operator = PotentialOperator::from(operator).unwrap(); +// let dtype = DType::from(dtype).unwrap(); +// match dtype { +// DType::C32 => helmholtz_potential_assembler_new_internal::( +// *(wavenumber as *const f32), +// operator, +// dtype, +// ), +// DType::C64 => helmholtz_potential_assembler_new_internal::( +// *(wavenumber as *const f64), +// operator, +// dtype, +// ), +// _ => { +// panic!("Invalid data type"); +// } +// } +// } +// } diff --git a/src/helmholtz.rs b/src/helmholtz.rs new file mode 100644 index 00000000..e83bdfa3 --- /dev/null +++ b/src/helmholtz.rs @@ -0,0 +1,146 @@ +//! Helmholtz operators + +pub mod assembler { + use green_kernels::{helmholtz_3d::Helmholtz3dKernel, types::GreenKernelEvalType}; + use rlst::{rlst_dynamic_array2, DynamicArray, MatrixInverse, RlstScalar}; + + use crate::{ + assembly::{ + boundary::{ + integrands::{ + AdjointDoubleLayerBoundaryIntegrand, BoundaryIntegrandScalarProduct, + BoundaryIntegrandSum, DoubleLayerBoundaryIntegrand, + HypersingularCurlCurlBoundaryIntegrand, + HypersingularNormalNormalBoundaryIntegrand, SingleLayerBoundaryIntegrand, + }, + BoundaryAssembler, BoundaryAssemblerOptions, + }, + kernels::KernelEvaluator, + }, + traits::{BoundaryAssembly, FunctionSpace}, + }; + + /// Assembler for the Helmholtz single layer operator. + pub fn helmholtz_single_layer< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new( + Helmholtz3dKernel::new(wavenumber), + GreenKernelEvalType::Value, + ); + + let assembler = + BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } + + /// Assembler for the Helmholtz double layer operator. + pub fn helmholtz_double_layer< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new( + Helmholtz3dKernel::new(wavenumber), + GreenKernelEvalType::ValueDeriv, + ); + + let assembler = + BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } + + /// Assembler for the Helmholtz adjoint double layer operator. + pub fn helmholtz_adjoint_double_layer< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new( + Helmholtz3dKernel::new(wavenumber), + GreenKernelEvalType::ValueDeriv, + ); + + let assembler = BoundaryAssembler::new( + AdjointDoubleLayerBoundaryIntegrand::new(), + kernel, + options, + 4, + 0, + ); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } + + /// Assembler for the Helmholtz hypersingular operator. + pub fn helmholtz_hypersingular< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new( + Helmholtz3dKernel::new(wavenumber), + GreenKernelEvalType::ValueDeriv, + ); + + let integrand = BoundaryIntegrandSum::new( + HypersingularCurlCurlBoundaryIntegrand::new(), + BoundaryIntegrandScalarProduct::new( + num::cast::(-wavenumber.powi(2)).unwrap(), + HypersingularNormalNormalBoundaryIntegrand::new(), + ), + ); + + let assembler = BoundaryAssembler::new(integrand, kernel, options, 4, 1); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } +} diff --git a/src/laplace.rs b/src/laplace.rs new file mode 100644 index 00000000..63d0c818 --- /dev/null +++ b/src/laplace.rs @@ -0,0 +1,126 @@ +//! Laplace operators + +pub mod assembler { + use green_kernels::{laplace_3d::Laplace3dKernel, types::GreenKernelEvalType}; + use rlst::{rlst_dynamic_array2, DynamicArray, MatrixInverse, RlstScalar}; + + use crate::{ + assembly::{ + boundary::{ + integrands::{ + AdjointDoubleLayerBoundaryIntegrand, DoubleLayerBoundaryIntegrand, + HypersingularCurlCurlBoundaryIntegrand, SingleLayerBoundaryIntegrand, + }, + BoundaryAssembler, BoundaryAssemblerOptions, + }, + kernels::KernelEvaluator, + }, + traits::{BoundaryAssembly, FunctionSpace}, + }; + + /// Assembler for the Laplace single layer operator. + pub fn laplace_single_layer< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value); + + let assembler = + BoundaryAssembler::new(SingleLayerBoundaryIntegrand::new(), kernel, options, 1, 0); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } + + /// Assembler for the Laplace double layer operator. + pub fn laplace_double_layer< + T: RlstScalar + MatrixInverse, + Space: FunctionSpace + Sync, + >( + 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]); + + let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); + + let assembler = + BoundaryAssembler::new(DoubleLayerBoundaryIntegrand::new(), kernel, options, 4, 0); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } + + /// 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, + options: BoundaryAssemblerOptions, + ) -> DynamicArray { + let nrows = trial_space.global_size(); + let ncols = test_space.global_size(); + + let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); + + let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); + + let assembler = BoundaryAssembler::new( + AdjointDoubleLayerBoundaryIntegrand::new(), + kernel, + options, + 4, + 0, + ); + + assembler.assemble_into_dense(&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, + options: BoundaryAssemblerOptions, + ) -> DynamicArray { + let nrows = trial_space.global_size(); + let ncols = test_space.global_size(); + + let mut output = rlst_dynamic_array2!(T, [nrows, ncols]); + + let kernel = KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::ValueDeriv); + + let assembler = BoundaryAssembler::new( + HypersingularCurlCurlBoundaryIntegrand::new(), + kernel, + options, + 4, + 1, + ); + + assembler.assemble_into_dense(&mut output, trial_space, test_space); + + output + } +} diff --git a/src/lib.rs b/src/lib.rs index a3850a08..1f327ba2 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -5,6 +5,8 @@ pub mod assembly; pub mod bindings; pub mod function; +pub mod helmholtz; +pub mod laplace; pub mod traits; #[cfg(test)] diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs index 51243c7c..4d2d9d44 100644 --- a/tests/compare_to_bempp_cl.rs +++ b/tests/compare_to_bempp_cl.rs @@ -1,7 +1,15 @@ use approx::*; -use bempp::assembly::{boundary::BoundaryAssembler, potential::PotentialAssembler}; +use bempp::assembly::boundary::BoundaryAssemblerOptions; +use bempp::assembly::potential::PotentialAssembler; use bempp::function::SerialFunctionSpace; -use bempp::traits::{BoundaryAssembly, FunctionSpace, PotentialAssembly}; +use bempp::helmholtz::assembler::{ + helmholtz_adjoint_double_layer, helmholtz_double_layer, helmholtz_hypersingular, + helmholtz_single_layer, +}; +use bempp::laplace::assembler::{ + laplace_adjoint_double_layer, laplace_double_layer, laplace_hypersingular, laplace_single_layer, +}; +use bempp::traits::{FunctionSpace, PotentialAssembly}; use cauchy::c64; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::Continuity; @@ -14,12 +22,7 @@ fn test_laplace_single_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - - let a = BoundaryAssembler::::new_laplace_single_layer(); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = laplace_single_layer(&space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -38,11 +41,7 @@ fn test_laplace_double_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = BoundaryAssembler::::new_laplace_double_layer(); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = laplace_double_layer(&space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -60,11 +59,7 @@ fn test_laplace_adjoint_double_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = BoundaryAssembler::::new_laplace_adjoint_double_layer(); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = laplace_adjoint_double_layer(&space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -77,36 +72,13 @@ fn test_laplace_adjoint_double_layer_dp0_dp0() { } } -#[test] -fn test_laplace_hypersingular_dp0_dp0() { - let grid = regular_sphere(0); - let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); - let space = SerialFunctionSpace::new(&grid, &element); - - let ndofs = space.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = BoundaryAssembler::::new_laplace_hypersingular(); - a.assemble_into_dense(&mut matrix, &space, &space); - - for i in 0..ndofs { - for j in 0..ndofs { - assert_relative_eq!(*matrix.get([i, j]).unwrap(), 0.0, epsilon = 1e-4); - } - } -} - #[test] fn test_laplace_hypersingular_p1_p1() { let grid = regular_sphere(0); let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = BoundaryAssembler::::new_laplace_hypersingular(); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = laplace_hypersingular(&space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -131,11 +103,7 @@ fn test_helmholtz_single_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]); - - let a = BoundaryAssembler::::new_helmholtz_single_layer(3.0); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = helmholtz_single_layer(3.0, &space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -154,11 +122,7 @@ fn test_helmholtz_double_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]); - - let a = BoundaryAssembler::::new_helmholtz_double_layer(3.0); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = helmholtz_double_layer(3.0, &space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -176,11 +140,8 @@ fn test_helmholtz_adjoint_double_layer_dp0_dp0() { let element = LagrangeElementFamily::::new(0, Continuity::Discontinuous); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]); - - let a = BoundaryAssembler::::new_helmholtz_adjoint_double_layer(3.0); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = + helmholtz_adjoint_double_layer(3.0, &space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip] @@ -199,11 +160,7 @@ fn test_helmholtz_hypersingular_p1_p1() { let element = LagrangeElementFamily::::new(1, Continuity::Standard); let space = SerialFunctionSpace::new(&grid, &element); - let ndofs = space.global_size(); - let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]); - - let a = BoundaryAssembler::::new_helmholtz_hypersingular(3.0); - a.assemble_into_dense(&mut matrix, &space, &space); + let matrix = helmholtz_hypersingular(3.0, &space, &space, BoundaryAssemblerOptions::default()); // Compare to result from bempp-cl #[rustfmt::skip]