From 83c792a3b94c60f54ad41a47214a2245ae2241d3 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 20 Aug 2024 17:09:11 +0100 Subject: [PATCH 01/13] working on refactor --- src/assembly.rs | 1 + src/assembly/boundary.rs | 252 +++++++++++++----- src/assembly/boundary/adjoint_double_layer.rs | 36 ++- src/assembly/boundary/double_layer.rs | 36 ++- src/assembly/boundary/hypersingular.rs | 88 ++++-- src/assembly/boundary/integrands.rs | 10 + .../integrands/adjoint_double_layer.rs | 75 ++++++ .../boundary/integrands/double_layer.rs | 85 ++++++ .../boundary/integrands/hypersingular.rs | 49 ++++ .../boundary/integrands/single_layer.rs | 53 ++++ src/assembly/boundary/single_layer.rs | 34 ++- src/assembly/common.rs | 40 +++ src/assembly/kernels.rs | 49 ++++ src/traits.rs | 2 + src/traits/assembly.rs | 76 ++++++ 15 files changed, 769 insertions(+), 117 deletions(-) create mode 100644 src/assembly/boundary/integrands.rs create mode 100644 src/assembly/boundary/integrands/adjoint_double_layer.rs create mode 100644 src/assembly/boundary/integrands/double_layer.rs create mode 100644 src/assembly/boundary/integrands/hypersingular.rs create mode 100644 src/assembly/boundary/integrands/single_layer.rs create mode 100644 src/assembly/kernels.rs create mode 100644 src/traits/assembly.rs diff --git a/src/assembly.rs b/src/assembly.rs index 7398c597..8a359bbb 100644 --- a/src/assembly.rs +++ b/src/assembly.rs @@ -2,6 +2,7 @@ pub mod boundary; pub(crate) mod common; pub mod fmm_tools; +pub mod kernels; pub mod potential; #[cfg(test)] diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index 48afc285..c2d9e081 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -2,6 +2,7 @@ pub(crate) mod adjoint_double_layer; pub(crate) mod double_layer; pub(crate) mod hypersingular; +mod integrands; pub(crate) mod single_layer; pub use adjoint_double_layer::AdjointDoubleLayerAssembler; @@ -9,15 +10,18 @@ pub use double_layer::DoubleLayerAssembler; pub use hypersingular::HypersingularAssembler; pub use single_layer::SingleLayerAssembler; -use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData}; +use crate::assembly::common::{ + equal_grids, AssemblerGeometry, RawData2D, RlstArray, SparseMatrixData, +}; use crate::quadrature::duffy::{ quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy, }; use crate::quadrature::simplex_rules::simplex_rule; use crate::quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; -use crate::traits::FunctionSpace; #[cfg(feature = "mpi")] use crate::traits::ParallelFunctionSpace; +use crate::traits::{BoundaryIntegrand, FunctionSpace, KernelEvaluator}; +use itertools::izip; #[cfg(feature = "mpi")] use mpi::traits::Communicator; use ndelement::reference_cell; @@ -25,10 +29,12 @@ use ndelement::traits::FiniteElement; use ndelement::types::ReferenceCellType; use ndgrid::traits::{Entity, GeometryMap, Grid, Topology}; use ndgrid::types::Ownership; +use num::Zero; use rayon::prelude::*; use rlst::{ - rlst_dynamic_array2, rlst_dynamic_array3, rlst_dynamic_array4, CsrMatrix, MatrixInverse, - RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, UnsafeRandomAccessByRef, + rlst_dynamic_array2, rlst_dynamic_array3, rlst_dynamic_array4, CsrMatrix, DefaultIterator, + DefaultIteratorMut, MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, + UnsafeRandomAccessByRef, }; use std::collections::HashMap; @@ -101,6 +107,148 @@ fn get_singular_quadrature_rule( } } +struct SingularCellPairAssembler< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, +> { + integrand: &'a I, + kernel: &'a K, + test_evaluator: G, + trial_evaluator: G, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + k: RlstArray, + test_mapped_pts: RlstArray, + trial_mapped_pts: RlstArray, + test_normals: RlstArray, + trial_normals: RlstArray, + test_jacobians: RlstArray, + trial_jacobians: RlstArray, + test_jdet: Vec, + trial_jdet: Vec, +} + +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, + > SingularCellPairAssembler<'a, T, I, G, K> +{ + pub fn new( + npts: usize, + deriv_size: usize, + integrand: &'a I, + kernel: &'a K, + test_evaluator: G, + trial_evaluator: G, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + ) -> Self { + Self { + integrand, + kernel, + test_evaluator, + trial_evaluator, + test_table, + trial_table, + k: rlst_dynamic_array2!(T, [deriv_size, npts]), + test_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), + trial_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), + test_normals: rlst_dynamic_array2!(T::Real, [3, npts]), + trial_normals: rlst_dynamic_array2!(T::Real, [3, npts]), + test_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), + trial_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), + test_jdet: vec![T::Real::zero(); npts], + trial_jdet: vec![T::Real::zero(); npts], + } + } +} +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, + > SingularCellPairAssembler<'a, T, I, G, K> +{ + // TODO make this a trait + fn assemble( + &mut self, + local_mat: &mut RlstArray, + weights: &[T::Real], + test_cell: usize, + trial_cell: usize, + ) { + self.test_evaluator + .points(test_cell, self.test_mapped_pts.data_mut()); + self.test_evaluator.jacobians_dets_normals( + test_cell, + self.test_jacobians.data_mut(), + &mut self.test_jdet, + self.test_normals.data_mut(), + ); + + self.trial_evaluator + .points(trial_cell, self.trial_mapped_pts.data_mut()); + self.trial_evaluator.jacobians_dets_normals( + trial_cell, + self.trial_jacobians.data_mut(), + &mut self.trial_jdet, + self.trial_normals.data_mut(), + ); + + // TODO: + /* + self.kernel.assemble_pairwise_st( + self.test_mapped_pts.data(), + self.trial_mapped_pts.data(), + self.k.data_mut(), + ); + */ + + let test_geometry = AssemblerGeometry::new( + &self.test_mapped_pts, + &self.test_normals, + &self.test_jacobians, + &self.test_jdet, + ); + let trial_geometry = AssemblerGeometry::new( + &self.trial_mapped_pts, + &self.trial_normals, + &self.trial_jacobians, + &self.trial_jdet, + ); + + for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { + for (test_i, entry) in col.iter_mut().enumerate() { + *entry = T::zero(); + for (index, wt) in weights.iter().enumerate() { + unsafe { + *entry += self.integrand.evaluate_singular( + self.test_table, + self.trial_table, + index, + test_i, + trial_i, + &self.k, + &test_geometry, + &trial_geometry, + ) * num::cast::( + *wt * *self.test_jdet.get_unchecked(index) + * *self.trial_jdet.get_unchecked(index), + ) + .unwrap(); + } + } + } + } + } +} /// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells #[allow(clippy::too_many_arguments)] fn assemble_batch_singular< @@ -138,82 +286,42 @@ fn assemble_batch_singular< assert_eq!(grid.geometry_dim(), 3); assert_eq!(grid.topology_dim(), 2); - // Memory assignment to be moved elsewhere as passed into here mutable? - let mut k = rlst_dynamic_array2!(T, [deriv_size, npts]); - let zero = num::cast::(0.0).unwrap(); - let mut test_jdet = vec![zero; npts]; - let mut test_mapped_pts = rlst_dynamic_array2!(T::Real, [3, npts]); - let mut test_jacobians = rlst_dynamic_array2!(T::Real, [6, npts]); - let mut test_normals = rlst_dynamic_array2!(T::Real, [3, npts]); - - let mut trial_jdet = vec![zero; npts]; - let mut trial_mapped_pts = rlst_dynamic_array2!(T::Real, [3, npts]); - let mut trial_jacobians = rlst_dynamic_array2!(T::Real, [6, npts]); - let mut trial_normals = rlst_dynamic_array2!(T::Real, [3, npts]); - let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); - for (test_cell, trial_cell) in cell_pairs { - test_evaluator.points(*test_cell, test_mapped_pts.data_mut()); - test_evaluator.jacobians_dets_normals( - *test_cell, - test_jacobians.data_mut(), - &mut test_jdet, - test_normals.data_mut(), - ); - - trial_evaluator.points(*trial_cell, trial_mapped_pts.data_mut()); - trial_evaluator.jacobians_dets_normals( - *trial_cell, - trial_jacobians.data_mut(), - &mut trial_jdet, - trial_normals.data_mut(), - ); + let mut a = SingularCellPairAssembler::new( + npts, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + ); - assembler.kernel_assemble_pairwise_st( - test_mapped_pts.data(), - trial_mapped_pts.data(), - k.data_mut(), - ); + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); + for (test_cell, trial_cell) in cell_pairs { + a.assemble(&mut local_mat, weights, *test_cell, *trial_cell); let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - for (test_i, test_dof) in test_dofs.iter().enumerate() { - for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { - let mut sum = num::cast::(0.0).unwrap(); - for (index, wt) in weights.iter().enumerate() { - unsafe { - sum += assembler.singular_kernel_value( - &k, - &test_normals, - &trial_normals, - index, - ) * assembler.test_trial_product( - test_table, - trial_table, - &test_jacobians, - &trial_jacobians, - &test_jdet, - &trial_jdet, - index, - index, - test_i, - trial_i, - ) * num::cast::( - *wt * *test_jdet.get_unchecked(index) - * *trial_jdet.get_unchecked(index), - ) - .unwrap(); - } - } + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { output.rows.push(test_space.global_dof_index(*test_dof)); output.cols.push(trial_space.global_dof_index(*trial_dof)); - output.data.push(sum); + output.data.push(entry); } } } + output } @@ -536,11 +644,21 @@ pub trait BoundaryAssembler: Sync + Sized { /// Scalar type type T: RlstScalar + MatrixInverse; + /// Integrand type + type Integrand: BoundaryIntegrand; + /// Kernel type + type Kernel: KernelEvaluator; /// Number of derivatives const DERIV_SIZE: usize; /// Number of derivatives needed in basis function tables const TABLE_DERIVS: usize; + /// Get integrand + fn integrand(&self) -> &Self::Integrand; + + /// Get integrand + fn kernel(&self) -> &Self::Kernel; + /// Get assembler options fn options(&self) -> &BoundaryAssemblerOptions; diff --git a/src/assembly/boundary/adjoint_double_layer.rs b/src/assembly/boundary/adjoint_double_layer.rs index 36f3d000..6772c8ef 100644 --- a/src/assembly/boundary/adjoint_double_layer.rs +++ b/src/assembly/boundary/adjoint_double_layer.rs @@ -1,18 +1,25 @@ //! Adjoint double layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; -use crate::assembly::common::{GreenKernelEvalType, RlstArray}; +use crate::assembly::{ + boundary::integrands::AdjointDoubleLayerBoundaryIntegrand, + common::{GreenKernelEvalType, RlstArray}, + kernels::KernelEvaluator, +}; +use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; /// Assembler for a adjoint double layer operator pub struct AdjointDoubleLayerAssembler> { - kernel: K, + integrand: AdjointDoubleLayerBoundaryIntegrand, + kernel: KernelEvaluator, options: BoundaryAssemblerOptions, } impl> AdjointDoubleLayerAssembler { /// Create a new adjoint double layer assembler - pub fn new(kernel: K) -> Self { + pub fn new(kernel: KernelEvaluator) -> Self { Self { + integrand: AdjointDoubleLayerBoundaryIntegrand::new(), kernel, options: BoundaryAssemblerOptions::default(), } @@ -21,7 +28,9 @@ impl> AdjointDoubleLayerAssemble impl AdjointDoubleLayerAssembler> { /// Create a new Laplace adjoint double layer assembler pub fn new_laplace() -> Self { - Self::new(Laplace3dKernel::::new()) + Self::new(KernelEvaluator::new_laplace( + GreenKernelEvalType::ValueDeriv, + )) } } impl + MatrixInverse> @@ -29,7 +38,10 @@ impl + MatrixInverse> { /// Create a new Helmholtz adjoint double layer assembler pub fn new_helmholtz(wavenumber: T::Real) -> Self { - Self::new(Helmholtz3dKernel::::new(wavenumber)) + Self::new(KernelEvaluator::new_helmholtz( + wavenumber, + GreenKernelEvalType::ValueDeriv, + )) } } @@ -39,6 +51,14 @@ impl> BoundaryAssembler const DERIV_SIZE: usize = 4; const TABLE_DERIVS: usize = 0; type T = T; + type Integrand = AdjointDoubleLayerBoundaryIntegrand; + type Kernel = KernelEvaluator; + fn integrand(&self) -> &AdjointDoubleLayerBoundaryIntegrand { + &self.integrand + } + fn kernel(&self) -> &KernelEvaluator { + &self.kernel + } fn options(&self) -> &BoundaryAssemblerOptions { &self.options } @@ -80,11 +100,9 @@ impl> BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::ValueDeriv, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::ValueDeriv, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } } diff --git a/src/assembly/boundary/double_layer.rs b/src/assembly/boundary/double_layer.rs index 51083622..cca63400 100644 --- a/src/assembly/boundary/double_layer.rs +++ b/src/assembly/boundary/double_layer.rs @@ -1,18 +1,25 @@ //! Double layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; -use crate::assembly::common::{GreenKernelEvalType, RlstArray}; +use crate::assembly::{ + boundary::integrands::DoubleLayerBoundaryIntegrand, + common::{GreenKernelEvalType, RlstArray}, + kernels::KernelEvaluator, +}; +use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; /// Assembler for a double layer operator pub struct DoubleLayerAssembler> { - kernel: K, + integrand: DoubleLayerBoundaryIntegrand, + kernel: KernelEvaluator, options: BoundaryAssemblerOptions, } impl> DoubleLayerAssembler { /// Create a new double layer assembler - pub fn new(kernel: K) -> Self { + pub fn new(kernel: KernelEvaluator) -> Self { Self { + integrand: DoubleLayerBoundaryIntegrand::new(), kernel, options: BoundaryAssemblerOptions::default(), } @@ -21,13 +28,18 @@ impl> DoubleLayerAssembler impl DoubleLayerAssembler> { /// Create a new Laplace double layer assembler pub fn new_laplace() -> Self { - Self::new(Laplace3dKernel::::new()) + Self::new(KernelEvaluator::new_laplace( + GreenKernelEvalType::ValueDeriv, + )) } } impl + MatrixInverse> DoubleLayerAssembler> { /// Create a new Helmholtz double layer assembler pub fn new_helmholtz(wavenumber: T::Real) -> Self { - Self::new(Helmholtz3dKernel::::new(wavenumber)) + Self::new(KernelEvaluator::new_helmholtz( + wavenumber, + GreenKernelEvalType::ValueDeriv, + )) } } @@ -37,6 +49,14 @@ impl> BoundaryAssembler const DERIV_SIZE: usize = 4; const TABLE_DERIVS: usize = 0; type T = T; + type Integrand = DoubleLayerBoundaryIntegrand; + type Kernel = KernelEvaluator; + fn integrand(&self) -> &DoubleLayerBoundaryIntegrand { + &self.integrand + } + fn kernel(&self) -> &KernelEvaluator { + &self.kernel + } fn options(&self) -> &BoundaryAssemblerOptions { &self.options } @@ -78,11 +98,9 @@ impl> BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::ValueDeriv, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::ValueDeriv, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } } diff --git a/src/assembly/boundary/hypersingular.rs b/src/assembly/boundary/hypersingular.rs index fbf9ae2a..dfc0cf50 100644 --- a/src/assembly/boundary/hypersingular.rs +++ b/src/assembly/boundary/hypersingular.rs @@ -1,7 +1,12 @@ //! Hypersingular assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; -use crate::assembly::common::{equal_grids, GreenKernelEvalType, RlstArray, SparseMatrixData}; +use crate::assembly::{ + boundary::integrands::HypersingularBoundaryIntegrand, + common::{equal_grids, GreenKernelEvalType, RlstArray, SparseMatrixData}, + kernels::KernelEvaluator, +}; use crate::traits::FunctionSpace; +use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use ndelement::traits::FiniteElement; use ndelement::types::ReferenceCellType; @@ -65,12 +70,12 @@ unsafe fn hyp_test_trial_product( /// Assembler for a hypersingular operator pub struct HypersingularAssembler> { - kernel: K, + kernel: KernelEvaluator, options: BoundaryAssemblerOptions, } impl> HypersingularAssembler { /// Create a new hypersingular assembler - pub fn new(kernel: K) -> Self { + pub fn new(kernel: KernelEvaluator) -> Self { Self { kernel, options: BoundaryAssemblerOptions::default(), @@ -80,13 +85,18 @@ impl> HypersingularAssembler HypersingularAssembler> { /// Create a new Laplace hypersingular assembler pub fn new_laplace() -> Self { - Self::new(Laplace3dKernel::::new()) + Self::new(KernelEvaluator::new_laplace( + GreenKernelEvalType::ValueDeriv, + )) } } impl + MatrixInverse> HypersingularAssembler> { /// Create a new Helmholtz hypersingular assembler pub fn new_helmholtz(wavenumber: T::Real) -> Self { - Self::new(Helmholtz3dKernel::::new(wavenumber)) + Self::new(KernelEvaluator::new_helmholtz( + wavenumber, + GreenKernelEvalType::ValueDeriv, + )) } } @@ -96,6 +106,14 @@ impl BoundaryAssembler const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 1; type T = T; + type Integrand = HypersingularBoundaryIntegrand; + type Kernel = KernelEvaluator>; + fn integrand(&self) -> &HypersingularBoundaryIntegrand { + panic!(); + } + fn kernel(&self) -> &KernelEvaluator> { + &self.kernel + } fn options(&self) -> &BoundaryAssemblerOptions { &self.options } @@ -127,12 +145,10 @@ impl BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } #[allow(clippy::too_many_arguments)] unsafe fn test_trial_product( @@ -165,14 +181,17 @@ impl BoundaryAssembler /// Assembler for curl-curl term of Helmholtz hypersingular operator struct HelmholtzHypersingularCurlCurlAssembler<'a, T: RlstScalar + MatrixInverse> { - kernel: &'a Helmholtz3dKernel, + kernel: &'a KernelEvaluator>, options: &'a BoundaryAssemblerOptions, } impl<'a, T: RlstScalar + MatrixInverse> HelmholtzHypersingularCurlCurlAssembler<'a, T> { /// Create a new assembler - pub fn new(kernel: &'a Helmholtz3dKernel, options: &'a BoundaryAssemblerOptions) -> Self { + pub fn new( + kernel: &'a KernelEvaluator>, + options: &'a BoundaryAssemblerOptions, + ) -> Self { Self { kernel, options } } } @@ -182,6 +201,14 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 1; type T = T; + type Integrand = HypersingularBoundaryIntegrand; + type Kernel = KernelEvaluator>; + fn integrand(&self) -> &HypersingularBoundaryIntegrand { + panic!(); + } + fn kernel(&self) -> &KernelEvaluator> { + panic!(); + } fn options(&self) -> &BoundaryAssemblerOptions { self.options } @@ -213,12 +240,10 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } #[allow(clippy::too_many_arguments)] unsafe fn test_trial_product( @@ -251,14 +276,17 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler /// Assembler for normal normal term of Helmholtz hypersingular boundary operator struct HelmholtzHypersingularNormalNormalAssembler<'a, T: RlstScalar + MatrixInverse> { - kernel: &'a Helmholtz3dKernel, + kernel: &'a KernelEvaluator>, options: &'a BoundaryAssemblerOptions, } impl<'a, T: RlstScalar + MatrixInverse> HelmholtzHypersingularNormalNormalAssembler<'a, T> { /// Create a new assembler - pub fn new(kernel: &'a Helmholtz3dKernel, options: &'a BoundaryAssemblerOptions) -> Self { + pub fn new( + kernel: &'a KernelEvaluator>, + options: &'a BoundaryAssemblerOptions, + ) -> Self { Self { kernel, options } } } @@ -268,6 +296,14 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 0; type T = T; + type Integrand = HypersingularBoundaryIntegrand; + type Kernel = KernelEvaluator>; + fn integrand(&self) -> &HypersingularBoundaryIntegrand { + panic!(); + } + fn kernel(&self) -> &KernelEvaluator> { + panic!(); + } fn options(&self) -> &BoundaryAssemblerOptions { self.options } @@ -283,7 +319,7 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler ) -> T { -*k.get_unchecked([0, index]) * num::cast::( - self.kernel.wavenumber.powi(2) + self.kernel.kernel.wavenumber.powi(2) * (*trial_normals.get_unchecked([0, index]) * *test_normals.get_unchecked([0, index]) + *trial_normals.get_unchecked([1, index]) @@ -301,7 +337,7 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler test_index: usize, trial_index: usize, ) -> T { - -num::cast::(self.kernel.wavenumber.powi(2)).unwrap() + -num::cast::(self.kernel.kernel.wavenumber.powi(2)).unwrap() * *k.get_unchecked([0, test_index, trial_index]) * (num::cast::(*trial_normals.get_unchecked([0, trial_index])).unwrap() * num::cast::(*test_normals.get_unchecked([0, test_index])).unwrap() @@ -318,12 +354,10 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } } @@ -333,6 +367,14 @@ impl + MatrixInverse> BoundaryAssembler const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 1; type T = T; + type Integrand = HypersingularBoundaryIntegrand; + type Kernel = KernelEvaluator>; + fn integrand(&self) -> &HypersingularBoundaryIntegrand { + panic!(); + } + fn kernel(&self) -> &KernelEvaluator> { + panic!(); + } fn options(&self) -> &BoundaryAssemblerOptions { &self.options } diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs new file mode 100644 index 00000000..c7751b6d --- /dev/null +++ b/src/assembly/boundary/integrands.rs @@ -0,0 +1,10 @@ +//! Integrands +mod adjoint_double_layer; +mod double_layer; +mod hypersingular; +mod single_layer; + +pub use adjoint_double_layer::AdjointDoubleLayerBoundaryIntegrand; +pub use double_layer::DoubleLayerBoundaryIntegrand; +pub use hypersingular::HypersingularBoundaryIntegrand; +pub use single_layer::SingleLayerBoundaryIntegrand; diff --git a/src/assembly/boundary/integrands/adjoint_double_layer.rs b/src/assembly/boundary/integrands/adjoint_double_layer.rs new file mode 100644 index 00000000..049dbfe1 --- /dev/null +++ b/src/assembly/boundary/integrands/adjoint_double_layer.rs @@ -0,0 +1,75 @@ +//! Adjoint double layer integrand +use crate::assembly::common::RlstArray; +use crate::traits::{BoundaryIntegrand, CellGeometry}; +use rlst::{RlstScalar, UnsafeRandomAccessByRef}; + +pub struct AdjointDoubleLayerBoundaryIntegrand { + _t: std::marker::PhantomData, +} + +impl AdjointDoubleLayerBoundaryIntegrand { + pub fn new() -> Self { + Self { + _t: std::marker::PhantomData, + } + } +} + +impl BoundaryIntegrand for AdjointDoubleLayerBoundaryIntegrand { + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + _trial_geometry: &impl CellGeometry, + ) -> T { + -(*k.get_unchecked([1, test_point_index, trial_point_index]) + * num::cast::( + *test_geometry.normals().get_unchecked([0, test_point_index]), + ) + .unwrap() + + *k.get_unchecked([2, test_point_index, trial_point_index]) + * num::cast::( + *test_geometry.normals().get_unchecked([1, test_point_index]), + ) + .unwrap() + + *k.get_unchecked([3, test_point_index, trial_point_index]) + * num::cast::( + *test_geometry.normals().get_unchecked([2, test_point_index]), + ) + .unwrap()) + * *test_table.get_unchecked([0, test_point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, trial_point_index, trial_basis_index, 0]) + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + _trial_geometry: &impl CellGeometry, + ) -> T { + -(*k.get_unchecked([1, point_index]) + * num::cast::(*test_geometry.normals().get_unchecked([0, point_index])) + .unwrap() + + *k.get_unchecked([2, point_index]) + * num::cast::(*test_geometry.normals().get_unchecked([1, point_index])) + .unwrap() + + *k.get_unchecked([3, point_index]) + * num::cast::(*test_geometry.normals().get_unchecked([2, point_index])) + .unwrap()) + * *test_table.get_unchecked([0, point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) + } +} diff --git a/src/assembly/boundary/integrands/double_layer.rs b/src/assembly/boundary/integrands/double_layer.rs new file mode 100644 index 00000000..835b497e --- /dev/null +++ b/src/assembly/boundary/integrands/double_layer.rs @@ -0,0 +1,85 @@ +//! Double layer integrand +use crate::assembly::common::RlstArray; +use crate::traits::{BoundaryIntegrand, CellGeometry}; +use rlst::{RlstScalar, UnsafeRandomAccessByRef}; + +pub struct DoubleLayerBoundaryIntegrand { + _t: std::marker::PhantomData, +} + +impl DoubleLayerBoundaryIntegrand { + pub fn new() -> Self { + Self { + _t: std::marker::PhantomData, + } + } +} + +impl BoundaryIntegrand for DoubleLayerBoundaryIntegrand { + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + (*k.get_unchecked([1, test_point_index, trial_point_index]) + * num::cast::( + *trial_geometry + .normals() + .get_unchecked([0, trial_point_index]), + ) + .unwrap() + + *k.get_unchecked([2, test_point_index, trial_point_index]) + * num::cast::( + *trial_geometry + .normals() + .get_unchecked([1, trial_point_index]), + ) + .unwrap() + + *k.get_unchecked([3, test_point_index, trial_point_index]) + * num::cast::( + *trial_geometry + .normals() + .get_unchecked([2, trial_point_index]), + ) + .unwrap()) + * *test_table.get_unchecked([0, test_point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, trial_point_index, trial_basis_index, 0]) + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + (*k.get_unchecked([1, point_index]) + * num::cast::(*trial_geometry.normals().get_unchecked([0, point_index])) + .unwrap() + + *k.get_unchecked([2, point_index]) + * num::cast::( + *trial_geometry.normals().get_unchecked([1, point_index]), + ) + .unwrap() + + *k.get_unchecked([3, point_index]) + * num::cast::( + *trial_geometry.normals().get_unchecked([2, point_index]), + ) + .unwrap()) + * *test_table.get_unchecked([0, point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) + } +} diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs new file mode 100644 index 00000000..38683829 --- /dev/null +++ b/src/assembly/boundary/integrands/hypersingular.rs @@ -0,0 +1,49 @@ +//! Hypersingular integrand +use crate::assembly::common::RlstArray; +use crate::traits::{BoundaryIntegrand, CellGeometry}; +use rlst::{RlstScalar, UnsafeRandomAccessByRef}; + +pub struct HypersingularBoundaryIntegrand { + _t: std::marker::PhantomData, +} + +impl HypersingularBoundaryIntegrand { + pub fn new() -> Self { + Self { + _t: std::marker::PhantomData, + } + } +} + +impl BoundaryIntegrand for HypersingularBoundaryIntegrand { + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + panic!(); + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + panic!(); + } +} diff --git a/src/assembly/boundary/integrands/single_layer.rs b/src/assembly/boundary/integrands/single_layer.rs new file mode 100644 index 00000000..2e9803c8 --- /dev/null +++ b/src/assembly/boundary/integrands/single_layer.rs @@ -0,0 +1,53 @@ +//! Single layer integrand +use crate::assembly::common::RlstArray; +use crate::traits::{BoundaryIntegrand, CellGeometry}; +use rlst::{RlstScalar, UnsafeRandomAccessByRef}; + +pub struct SingleLayerBoundaryIntegrand { + _t: std::marker::PhantomData, +} + +impl SingleLayerBoundaryIntegrand { + pub fn new() -> Self { + Self { + _t: std::marker::PhantomData, + } + } +} + +impl BoundaryIntegrand for SingleLayerBoundaryIntegrand { + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + _trial_geometry: &impl CellGeometry, + ) -> T { + *k.get_unchecked([0, test_point_index, trial_point_index]) + * *test_table.get_unchecked([0, test_point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, trial_point_index, trial_basis_index, 0]) + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + _test_geometry: &impl CellGeometry, + _trial_geometry: &impl CellGeometry, + ) -> T { + *k.get_unchecked([0, point_index]) + * *test_table.get_unchecked([0, point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) + } +} diff --git a/src/assembly/boundary/single_layer.rs b/src/assembly/boundary/single_layer.rs index e3f812c8..425d2ade 100644 --- a/src/assembly/boundary/single_layer.rs +++ b/src/assembly/boundary/single_layer.rs @@ -1,18 +1,25 @@ //! Single layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; -use crate::assembly::common::{GreenKernelEvalType, RlstArray}; +use crate::assembly::{ + boundary::integrands::SingleLayerBoundaryIntegrand, + common::{GreenKernelEvalType, RlstArray}, + kernels::KernelEvaluator, +}; +use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; /// Assembler for a single layer operator pub struct SingleLayerAssembler> { - kernel: K, + integrand: SingleLayerBoundaryIntegrand, + kernel: KernelEvaluator, options: BoundaryAssemblerOptions, } impl> SingleLayerAssembler { /// Create a new single layer assembler - pub fn new(kernel: K) -> Self { + pub fn new(kernel: KernelEvaluator) -> Self { Self { + integrand: SingleLayerBoundaryIntegrand::new(), kernel, options: BoundaryAssemblerOptions::default(), } @@ -21,13 +28,16 @@ impl> SingleLayerAssembler impl SingleLayerAssembler> { /// Create a new Laplace single layer assembler pub fn new_laplace() -> Self { - Self::new(Laplace3dKernel::::new()) + Self::new(KernelEvaluator::new_laplace(GreenKernelEvalType::Value)) } } impl + MatrixInverse> SingleLayerAssembler> { /// Create a new Helmholtz single layer assembler pub fn new_helmholtz(wavenumber: T::Real) -> Self { - Self::new(Helmholtz3dKernel::::new(wavenumber)) + Self::new(KernelEvaluator::new_helmholtz( + wavenumber, + GreenKernelEvalType::Value, + )) } } @@ -37,6 +47,14 @@ impl> BoundaryAssembler const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 0; type T = T; + type Integrand = SingleLayerBoundaryIntegrand; + type Kernel = KernelEvaluator; + fn integrand(&self) -> &SingleLayerBoundaryIntegrand { + &self.integrand + } + fn kernel(&self) -> &KernelEvaluator { + &self.kernel + } fn options(&self) -> &BoundaryAssemblerOptions { &self.options } @@ -68,11 +86,9 @@ impl> BoundaryAssembler targets: &[T::Real], result: &mut [T], ) { - self.kernel - .assemble_pairwise_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_pairwise_st(sources, targets, result); } fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel - .assemble_st(GreenKernelEvalType::Value, sources, targets, result); + self.kernel.assemble_st(sources, targets, result); } } diff --git a/src/assembly/common.rs b/src/assembly/common.rs index fb3db393..fe184e06 100644 --- a/src/assembly/common.rs +++ b/src/assembly/common.rs @@ -1,4 +1,5 @@ //! Common utility functions +use crate::traits::CellGeometry; pub(crate) use green_kernels::types::EvalType as GreenKernelEvalType; use ndgrid::traits::Grid; use rlst::{Array, BaseArray, MatrixInverse, RlstScalar, VectorContainer}; @@ -78,3 +79,42 @@ impl SparseMatrixData { } unsafe impl Sync for SparseMatrixData {} + +pub(crate) struct AssemblerGeometry<'a, T: RlstScalar> { + points: &'a RlstArray, + normals: &'a RlstArray, + jacobians: &'a RlstArray, + jdets: &'a [T], +} + +impl<'a, T: RlstScalar> AssemblerGeometry<'a, T> { + pub(crate) fn new( + points: &'a RlstArray, + normals: &'a RlstArray, + jacobians: &'a RlstArray, + jdets: &'a [T], + ) -> Self { + Self { + points, + normals, + jacobians, + jdets, + } + } +} + +impl<'a, T: RlstScalar> CellGeometry for AssemblerGeometry<'a, T> { + type T = T; + fn points(&self) -> &RlstArray { + self.points + } + fn normals(&self) -> &RlstArray { + self.normals + } + fn jacobians(&self) -> &RlstArray { + self.jacobians + } + fn jdets(&self) -> &[T] { + self.jdets + } +} diff --git a/src/assembly/kernels.rs b/src/assembly/kernels.rs new file mode 100644 index 00000000..a4337d5f --- /dev/null +++ b/src/assembly/kernels.rs @@ -0,0 +1,49 @@ +//! Green's function kernels +use crate::assembly::common::GreenKernelEvalType; +use crate::traits::KernelEvaluator as KernelEvaluatorTrait; +use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; +use rlst::RlstScalar; + +pub struct KernelEvaluator> { + pub(crate) kernel: K, + eval_type: GreenKernelEvalType, +} + +impl> KernelEvaluator { + pub fn new(kernel: K, eval_type: GreenKernelEvalType) -> Self { + Self { kernel, eval_type } + } +} +impl KernelEvaluator> { + pub fn new_laplace(eval_type: GreenKernelEvalType) -> Self { + Self::new(Laplace3dKernel::::new(), eval_type) + } +} +impl> KernelEvaluator> { + pub fn new_helmholtz(wavenumber: T::Real, eval_type: GreenKernelEvalType) -> Self { + Self::new(Helmholtz3dKernel::::new(wavenumber), eval_type) + } +} +impl> KernelEvaluatorTrait for KernelEvaluator { + type T = T; + + fn assemble_pairwise_st( + &self, + sources: &[::Real], + targets: &[::Real], + result: &mut [Self::T], + ) { + self.kernel + .assemble_pairwise_st(self.eval_type, sources, targets, result); + } + + fn assemble_st( + &self, + sources: &[::Real], + targets: &[::Real], + result: &mut [Self::T], + ) { + self.kernel + .assemble_st(self.eval_type, sources, targets, result); + } +} diff --git a/src/traits.rs b/src/traits.rs index 6930b5c0..30a4ec63 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -1,7 +1,9 @@ //! Trait definitions +mod assembly; mod function; +pub use assembly::{BoundaryIntegrand, CellGeometry, KernelEvaluator}; pub use function::FunctionSpace; #[cfg(feature = "mpi")] pub use function::ParallelFunctionSpace; diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs new file mode 100644 index 00000000..5fa6440e --- /dev/null +++ b/src/traits/assembly.rs @@ -0,0 +1,76 @@ +//! Assembly +use crate::assembly::common::RlstArray; +use rlst::RlstScalar; + +pub trait CellGeometry { + //! Cell geometry + /// Scalar type + type T: RlstScalar; + /// Points + fn points(&self) -> &RlstArray; + /// Normals + fn normals(&self) -> &RlstArray; + /// Jacobians + fn jacobians(&self) -> &RlstArray; + /// Determinants of jacobians + fn jdets(&self) -> &[Self::T]; +} + +pub trait BoundaryIntegrand { + //! Integrand + /// Scalar type + type T: RlstScalar; + + /// Evaluate integrand for a singular quadrature rule + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry::Real>, + trial_geometry: &impl CellGeometry::Real>, + ) -> Self::T; + + /// Evaluate integrand for a non-singular quadrature rule + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry::Real>, + trial_geometry: &impl CellGeometry::Real>, + ) -> Self::T; +} + +pub trait KernelEvaluator { + //! Kernel evaluator + /// Scalar type + type T: RlstScalar; + + /// Evaluate the kernel values for all source and target pairs + /// + /// For each source, the kernel is evaluated for exactly one target. This is equivalent to taking the diagonal of the matrix assembled by `assemble_st` + fn assemble_pairwise_st( + &self, + sources: &[::Real], + targets: &[::Real], + result: &mut [Self::T], + ); + + /// Evaluate the kernel values for all sources and all targets + /// + /// For every source, the kernel is evaluated for every target. + fn assemble_st( + &self, + sources: &[::Real], + targets: &[::Real], + result: &mut [Self::T], + ); +} From 144ada387fa5ef16f5e6725aa07a40638a34699f Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Tue, 20 Aug 2024 17:28:07 +0100 Subject: [PATCH 02/13] remove some no-longer-needed functions --- src/assembly/boundary.rs | 145 +++----- src/assembly/boundary/adjoint_double_layer.rs | 46 +-- src/assembly/boundary/double_layer.rs | 46 +-- src/assembly/boundary/hypersingular.rs | 317 +----------------- .../boundary/integrands/hypersingular.rs | 2 +- src/assembly/boundary/single_layer.rs | 36 +- 6 files changed, 51 insertions(+), 541 deletions(-) diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index c2d9e081..9e8f84fe 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -34,7 +34,6 @@ use rayon::prelude::*; use rlst::{ rlst_dynamic_array2, rlst_dynamic_array3, rlst_dynamic_array4, CsrMatrix, DefaultIterator, DefaultIteratorMut, MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, - UnsafeRandomAccessByRef, }; use std::collections::HashMap; @@ -202,14 +201,11 @@ impl< self.trial_normals.data_mut(), ); - // TODO: - /* self.kernel.assemble_pairwise_st( self.test_mapped_pts.data(), self.trial_mapped_pts.data(), self.k.data_mut(), ); - */ let test_geometry = AssemblerGeometry::new( &self.test_mapped_pts, @@ -409,7 +405,7 @@ fn assemble_batch_nonadjacent< continue; } - assembler.kernel_assemble_st( + assembler.kernel().assemble_st( test_mapped_pts.data(), trial_mapped_pts[trial_cell_i].data(), k.data_mut(), @@ -418,6 +414,19 @@ fn assemble_batch_nonadjacent< let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); + let test_geometry = AssemblerGeometry::new( + &test_mapped_pts, + &test_normals, + &test_jacobians, + &test_jdet, + ); + let trial_geometry = AssemblerGeometry::new( + &trial_mapped_pts[trial_cell_i], + &trial_normals[trial_cell_i], + &trial_jacobians[trial_cell_i], + &trial_jdet[trial_cell_i], + ); + for (test_i, test_dof) in test_dofs.iter().enumerate() { for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { for (trial_index, trial_wt) in trial_weights.iter().enumerate() { @@ -432,26 +441,18 @@ fn assemble_batch_nonadjacent< num::cast::(*test_wt * test_jdet[test_index]).unwrap(); for trial_index in 0..npts_trial { sum += unsafe { - assembler.nonsingular_kernel_value( - &k, - &test_normals, - &trial_normals[trial_cell_i], + assembler.integrand().evaluate_nonsingular( + test_table, + trial_table, test_index, trial_index, + test_i, + trial_i, + &k, + &test_geometry, + &trial_geometry, ) * test_integrand * *trial_integrands.get_unchecked(trial_index) - * assembler.test_trial_product( - test_table, - trial_table, - &test_jacobians, - &trial_jacobians[trial_cell_i], - &test_jdet, - &trial_jdet[trial_cell_i], - test_index, - trial_index, - test_i, - trial_i, - ) }; } } @@ -540,12 +541,21 @@ fn assemble_batch_singular_correction< trial_normals.data_mut(), ); - assembler.kernel_assemble_st( + assembler.kernel().assemble_st( test_mapped_pts.data(), trial_mapped_pts.data(), k.data_mut(), ); + let test_geometry = + AssemblerGeometry::new(&test_mapped_pts, &test_normals, &test_jacobians, &test_jdet); + let trial_geometry = AssemblerGeometry::new( + &trial_mapped_pts, + &trial_normals, + &trial_jacobians, + &trial_jdet, + ); + let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); for (test_i, test_dof) in test_dofs.iter().enumerate() { @@ -560,26 +570,18 @@ fn assemble_batch_singular_correction< num::cast::(*test_wt * test_jdet[test_index]).unwrap(); for trial_index in 0..npts_trial { sum += unsafe { - assembler.nonsingular_kernel_value( - &k, - &test_normals, - &trial_normals, + assembler.integrand().evaluate_nonsingular( + test_table, + trial_table, test_index, trial_index, + test_i, + trial_i, + &k, + &test_geometry, + &trial_geometry, ) * test_integrand * *trial_integrands.get_unchecked(trial_index) - * assembler.test_trial_product( - test_table, - trial_table, - &test_jacobians, - &trial_jacobians, - &test_jdet, - &trial_jdet, - test_index, - trial_index, - test_i, - trial_i, - ) }; } } @@ -692,73 +694,6 @@ pub trait BoundaryAssembler: Sync + Sized { self.options_mut().batch_size = size; } - /// Return the kernel value to use in the integrand when using a singular quadrature rule - /// - /// # Safety - /// This method is unsafe to allow `get_unchecked` to be used - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray<::Real, 2>, - trial_normals: &RlstArray<::Real, 2>, - index: usize, - ) -> Self::T; - - /// Return the kernel value to use in the integrand when using a non-singular quadrature rule - /// - /// # Safety - /// This method is unsafe to allow `get_unchecked` to be used - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray<::Real, 2>, - trial_normals: &RlstArray<::Real, 2>, - test_index: usize, - trial_index: usize, - ) -> Self::T; - - /// Evaluate the kernel values for all source and target pairs - /// - /// For each source, the kernel is evaluated for exactly one target. This is equivalent to taking the diagonal of the matrix assembled by `kernel_assemble_st` - fn kernel_assemble_pairwise_st( - &self, - sources: &[::Real], - targets: &[::Real], - result: &mut [Self::T], - ); - - /// Evaluate the kernel values for all sources and all targets - /// - /// For every source, the kernel is evaluated for every target. - fn kernel_assemble_st( - &self, - sources: &[::Real], - targets: &[::Real], - result: &mut [Self::T], - ); - - /// The product of a test and trial function - /// - /// # Safety - /// This method is unsafe to allow `get_unchecked` to be used - #[allow(clippy::too_many_arguments)] - unsafe fn test_trial_product( - &self, - test_table: &RlstArray, - trial_table: &RlstArray, - _test_jacobians: &RlstArray<::Real, 2>, - _trial_jacobians: &RlstArray<::Real, 2>, - _test_jdets: &[::Real], - _trial_jdets: &[::Real], - test_point_index: usize, - trial_point_index: usize, - test_basis_index: usize, - trial_basis_index: usize, - ) -> Self::T { - *test_table.get_unchecked([0, test_point_index, test_basis_index, 0]) - * *trial_table.get_unchecked([0, trial_point_index, trial_basis_index, 0]) - } - /// Assemble the singular contributions fn assemble_singular< TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, diff --git a/src/assembly/boundary/adjoint_double_layer.rs b/src/assembly/boundary/adjoint_double_layer.rs index 6772c8ef..ee6ed46b 100644 --- a/src/assembly/boundary/adjoint_double_layer.rs +++ b/src/assembly/boundary/adjoint_double_layer.rs @@ -1,13 +1,11 @@ //! Adjoint double layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; use crate::assembly::{ - boundary::integrands::AdjointDoubleLayerBoundaryIntegrand, - common::{GreenKernelEvalType, RlstArray}, + boundary::integrands::AdjointDoubleLayerBoundaryIntegrand, common::GreenKernelEvalType, kernels::KernelEvaluator, }; -use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; +use rlst::{MatrixInverse, RlstScalar}; /// Assembler for a adjoint double layer operator pub struct AdjointDoubleLayerAssembler> { @@ -65,44 +63,4 @@ impl> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { &mut self.options } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray, - _trial_normals: &RlstArray, - index: usize, - ) -> T { - -*k.get_unchecked([1, index]) - * num::cast::(*test_normals.get_unchecked([0, index])).unwrap() - - *k.get_unchecked([2, index]) - * num::cast::(*test_normals.get_unchecked([1, index])).unwrap() - - *k.get_unchecked([3, index]) - * num::cast::(*test_normals.get_unchecked([2, index])).unwrap() - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray, - _trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - -*k.get_unchecked([1, test_index, trial_index]) - * num::cast::(*test_normals.get_unchecked([0, test_index])).unwrap() - - *k.get_unchecked([2, test_index, trial_index]) - * num::cast::(*test_normals.get_unchecked([1, test_index])).unwrap() - - *k.get_unchecked([3, test_index, trial_index]) - * num::cast::(*test_normals.get_unchecked([2, test_index])).unwrap() - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } } diff --git a/src/assembly/boundary/double_layer.rs b/src/assembly/boundary/double_layer.rs index cca63400..8e441cd2 100644 --- a/src/assembly/boundary/double_layer.rs +++ b/src/assembly/boundary/double_layer.rs @@ -1,13 +1,11 @@ //! Double layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; use crate::assembly::{ - boundary::integrands::DoubleLayerBoundaryIntegrand, - common::{GreenKernelEvalType, RlstArray}, + boundary::integrands::DoubleLayerBoundaryIntegrand, common::GreenKernelEvalType, kernels::KernelEvaluator, }; -use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; +use rlst::{MatrixInverse, RlstScalar}; /// Assembler for a double layer operator pub struct DoubleLayerAssembler> { @@ -63,44 +61,4 @@ impl> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { &mut self.options } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - trial_normals: &RlstArray, - index: usize, - ) -> T { - *k.get_unchecked([1, index]) - * num::cast::(*trial_normals.get_unchecked([0, index])).unwrap() - + *k.get_unchecked([2, index]) - * num::cast::(*trial_normals.get_unchecked([1, index])).unwrap() - + *k.get_unchecked([3, index]) - * num::cast::(*trial_normals.get_unchecked([2, index])).unwrap() - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - *k.get_unchecked([1, test_index, trial_index]) - * num::cast::(*trial_normals.get_unchecked([0, trial_index])).unwrap() - + *k.get_unchecked([2, test_index, trial_index]) - * num::cast::(*trial_normals.get_unchecked([1, trial_index])).unwrap() - + *k.get_unchecked([3, test_index, trial_index]) - * num::cast::(*trial_normals.get_unchecked([2, trial_index])).unwrap() - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } } diff --git a/src/assembly/boundary/hypersingular.rs b/src/assembly/boundary/hypersingular.rs index dfc0cf50..0650e2d8 100644 --- a/src/assembly/boundary/hypersingular.rs +++ b/src/assembly/boundary/hypersingular.rs @@ -1,19 +1,13 @@ //! Hypersingular assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; use crate::assembly::{ - boundary::integrands::HypersingularBoundaryIntegrand, - common::{equal_grids, GreenKernelEvalType, RlstArray, SparseMatrixData}, + boundary::integrands::HypersingularBoundaryIntegrand, common::GreenKernelEvalType, kernels::KernelEvaluator, }; -use crate::traits::FunctionSpace; -use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use ndelement::traits::FiniteElement; -use ndelement::types::ReferenceCellType; -use ndgrid::traits::Grid; -use rlst::{MatrixInverse, RlstScalar, Shape, UnsafeRandomAccessByRef}; -use std::collections::HashMap; +use rlst::{MatrixInverse, RlstScalar}; +/* #[allow(clippy::too_many_arguments)] unsafe fn hyp_test_trial_product( test_table: &RlstArray, @@ -67,7 +61,7 @@ unsafe fn hyp_test_trial_product( / num::cast::(test_jdets[test_point_index] * trial_jdets[trial_point_index]) .unwrap() } - +*/ /// Assembler for a hypersingular operator pub struct HypersingularAssembler> { kernel: KernelEvaluator, @@ -120,63 +114,6 @@ impl BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { &mut self.options } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - index: usize, - ) -> T { - *k.get_unchecked([0, index]) - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - *k.get_unchecked([0, test_index, trial_index]) - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } - #[allow(clippy::too_many_arguments)] - unsafe fn test_trial_product( - &self, - test_table: &RlstArray, - trial_table: &RlstArray, - test_jacobians: &RlstArray, - trial_jacobians: &RlstArray, - test_jdets: &[T::Real], - trial_jdets: &[T::Real], - test_point_index: usize, - trial_point_index: usize, - test_basis_index: usize, - trial_basis_index: usize, - ) -> T { - hyp_test_trial_product::( - test_table, - trial_table, - test_jacobians, - trial_jacobians, - test_jdets, - trial_jdets, - test_point_index, - trial_point_index, - test_basis_index, - trial_basis_index, - ) - } } /// Assembler for curl-curl term of Helmholtz hypersingular operator @@ -215,63 +152,6 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { panic!("Cannot get mutable options") } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - index: usize, - ) -> T { - *k.get_unchecked([0, index]) - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - *k.get_unchecked([0, test_index, trial_index]) - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } - #[allow(clippy::too_many_arguments)] - unsafe fn test_trial_product( - &self, - test_table: &RlstArray, - trial_table: &RlstArray, - test_jacobians: &RlstArray, - trial_jacobians: &RlstArray, - test_jdets: &[T::Real], - trial_jdets: &[T::Real], - test_point_index: usize, - trial_point_index: usize, - test_basis_index: usize, - trial_basis_index: usize, - ) -> T { - hyp_test_trial_product::( - test_table, - trial_table, - test_jacobians, - trial_jacobians, - test_jdets, - trial_jdets, - test_point_index, - trial_point_index, - test_basis_index, - trial_basis_index, - ) - } } /// Assembler for normal normal term of Helmholtz hypersingular boundary operator @@ -310,55 +190,6 @@ impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { panic!("Cannot get mutable options") } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray, - trial_normals: &RlstArray, - index: usize, - ) -> T { - -*k.get_unchecked([0, index]) - * num::cast::( - self.kernel.kernel.wavenumber.powi(2) - * (*trial_normals.get_unchecked([0, index]) - * *test_normals.get_unchecked([0, index]) - + *trial_normals.get_unchecked([1, index]) - * *test_normals.get_unchecked([1, index]) - + *trial_normals.get_unchecked([2, index]) - * *test_normals.get_unchecked([2, index])), - ) - .unwrap() - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - test_normals: &RlstArray, - trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - -num::cast::(self.kernel.kernel.wavenumber.powi(2)).unwrap() - * *k.get_unchecked([0, test_index, trial_index]) - * (num::cast::(*trial_normals.get_unchecked([0, trial_index])).unwrap() - * num::cast::(*test_normals.get_unchecked([0, test_index])).unwrap() - + num::cast::(*trial_normals.get_unchecked([1, trial_index])).unwrap() - * num::cast::(*test_normals.get_unchecked([1, test_index])) - .unwrap() - + num::cast::(*trial_normals.get_unchecked([2, trial_index])).unwrap() - * num::cast::(*test_normals.get_unchecked([2, test_index])) - .unwrap()) - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } } impl + MatrixInverse> BoundaryAssembler @@ -381,144 +212,4 @@ impl + MatrixInverse> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { &mut self.options } - - unsafe fn singular_kernel_value( - &self, - _k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - _index: usize, - ) -> T { - panic!("Cannot directly use HypersingularAssembler for Helmholtz"); - } - unsafe fn nonsingular_kernel_value( - &self, - _k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - _test_index: usize, - _trial_index: usize, - ) -> T { - panic!("Cannot directly use HypersingularAssembler for Helmholtz"); - } - fn kernel_assemble_pairwise_st( - &self, - _sources: &[T::Real], - _targets: &[T::Real], - _result: &mut [T], - ) { - panic!("Cannot directly use HypersingularAssembler for Helmholtz"); - } - fn kernel_assemble_st(&self, _sources: &[T::Real], _targets: &[T::Real], _result: &mut [T]) { - panic!("Cannot directly use HypersingularAssembler for Helmholtz"); - } - - fn assemble_singular< - TestGrid: Grid + Sync, - TrialGrid: Grid + Sync, - Element: FiniteElement + Sync, - >( - &self, - shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> SparseMatrixData { - let curl_curl_assembler = - HelmholtzHypersingularCurlCurlAssembler::::new(&self.kernel, &self.options); - let normal_normal_assembler = - HelmholtzHypersingularNormalNormalAssembler::::new(&self.kernel, &self.options); - - let mut curlcurl = curl_curl_assembler.assemble_singular::( - shape, - trial_space, - test_space, - ); - curlcurl.add( - normal_normal_assembler.assemble_singular::( - shape, - trial_space, - test_space, - ), - ); - curlcurl - } - - fn assemble_singular_correction< - TestGrid: Grid + Sync, - TrialGrid: Grid + Sync, - Element: FiniteElement + Sync, - >( - &self, - shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> SparseMatrixData { - if !equal_grids(test_space.grid(), trial_space.grid()) { - // If the test and trial grids are different, there are no neighbouring triangles - return SparseMatrixData::new(shape); - } - - let curl_curl_assembler = - HelmholtzHypersingularCurlCurlAssembler::::new(&self.kernel, &self.options); - let normal_normal_assembler = - HelmholtzHypersingularNormalNormalAssembler::::new(&self.kernel, &self.options); - - let mut curlcurl = curl_curl_assembler - .assemble_singular_correction::( - shape, - trial_space, - test_space, - ); - curlcurl.add( - normal_normal_assembler.assemble_singular_correction::( - shape, - trial_space, - test_space, - ), - ); - curlcurl - } - - #[allow(clippy::too_many_arguments)] - fn assemble_nonsingular_into_dense< - TestGrid: Grid + Sync, - TrialGrid: Grid + Sync, - Element: FiniteElement + Sync, - >( - &self, - output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - trial_colouring: &HashMap>>, - test_colouring: &HashMap>>, - ) { - if !trial_space.is_serial() || !test_space.is_serial() { - panic!("Dense assembly can only be used for function spaces stored in serial"); - } - if output.shape()[0] != test_space.global_size() - || output.shape()[1] != trial_space.global_size() - { - panic!("Matrix has wrong shape"); - } - - let curl_curl_assembler = - HelmholtzHypersingularCurlCurlAssembler::::new(&self.kernel, &self.options); - let normal_normal_assembler = - HelmholtzHypersingularNormalNormalAssembler::::new(&self.kernel, &self.options); - - curl_curl_assembler.assemble_nonsingular_into_dense::( - output, - trial_space, - test_space, - trial_colouring, - test_colouring, - ); - normal_normal_assembler.assemble_nonsingular_into_dense::( - output, - trial_space, - test_space, - trial_colouring, - test_colouring, - ); - } } diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs index 38683829..933b976d 100644 --- a/src/assembly/boundary/integrands/hypersingular.rs +++ b/src/assembly/boundary/integrands/hypersingular.rs @@ -1,7 +1,7 @@ //! Hypersingular integrand use crate::assembly::common::RlstArray; use crate::traits::{BoundaryIntegrand, CellGeometry}; -use rlst::{RlstScalar, UnsafeRandomAccessByRef}; +use rlst::RlstScalar; pub struct HypersingularBoundaryIntegrand { _t: std::marker::PhantomData, diff --git a/src/assembly/boundary/single_layer.rs b/src/assembly/boundary/single_layer.rs index 425d2ade..4bf9f81a 100644 --- a/src/assembly/boundary/single_layer.rs +++ b/src/assembly/boundary/single_layer.rs @@ -1,13 +1,11 @@ //! Single layer assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; use crate::assembly::{ - boundary::integrands::SingleLayerBoundaryIntegrand, - common::{GreenKernelEvalType, RlstArray}, + boundary::integrands::SingleLayerBoundaryIntegrand, common::GreenKernelEvalType, kernels::KernelEvaluator, }; -use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; -use rlst::{MatrixInverse, RlstScalar, UnsafeRandomAccessByRef}; +use rlst::{MatrixInverse, RlstScalar}; /// Assembler for a single layer operator pub struct SingleLayerAssembler> { @@ -61,34 +59,4 @@ impl> BoundaryAssembler fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { &mut self.options } - unsafe fn singular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - index: usize, - ) -> T { - *k.get_unchecked([0, index]) - } - unsafe fn nonsingular_kernel_value( - &self, - k: &RlstArray, - _test_normals: &RlstArray, - _trial_normals: &RlstArray, - test_index: usize, - trial_index: usize, - ) -> T { - *k.get_unchecked([0, test_index, trial_index]) - } - fn kernel_assemble_pairwise_st( - &self, - sources: &[T::Real], - targets: &[T::Real], - result: &mut [T], - ) { - self.kernel.assemble_pairwise_st(sources, targets, result); - } - fn kernel_assemble_st(&self, sources: &[T::Real], targets: &[T::Real], result: &mut [T]) { - self.kernel.assemble_st(sources, targets, result); - } } From 9d292aa6d2fe9b8c193003a784e993a9083abe35 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 09:35:07 +0100 Subject: [PATCH 03/13] nonsingular cell pair assembler --- src/assembly/boundary.rs | 269 +++------------ src/assembly/boundary/cell_pair_assemblers.rs | 315 ++++++++++++++++++ src/assembly/kernels.rs | 4 + src/traits.rs | 2 +- src/traits/assembly.rs | 10 + 5 files changed, 368 insertions(+), 232 deletions(-) create mode 100644 src/assembly/boundary/cell_pair_assemblers.rs diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index 9e8f84fe..dcae5613 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -1,5 +1,6 @@ //! Assembly of boundary operators pub(crate) mod adjoint_double_layer; +mod cell_pair_assemblers; pub(crate) mod double_layer; pub(crate) mod hypersingular; mod integrands; @@ -20,7 +21,8 @@ use crate::quadrature::simplex_rules::simplex_rule; use crate::quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; #[cfg(feature = "mpi")] use crate::traits::ParallelFunctionSpace; -use crate::traits::{BoundaryIntegrand, FunctionSpace, KernelEvaluator}; +use crate::traits::{BoundaryIntegrand, CellPairAssembler, FunctionSpace, KernelEvaluator}; +use cell_pair_assemblers::{NonsingularCellPairAssembler, SingularCellPairAssembler}; use itertools::izip; #[cfg(feature = "mpi")] use mpi::traits::Communicator; @@ -29,11 +31,10 @@ use ndelement::traits::FiniteElement; use ndelement::types::ReferenceCellType; use ndgrid::traits::{Entity, GeometryMap, Grid, Topology}; use ndgrid::types::Ownership; -use num::Zero; use rayon::prelude::*; use rlst::{ rlst_dynamic_array2, rlst_dynamic_array3, rlst_dynamic_array4, CsrMatrix, DefaultIterator, - DefaultIteratorMut, MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, + MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, }; use std::collections::HashMap; @@ -106,145 +107,6 @@ fn get_singular_quadrature_rule( } } -struct SingularCellPairAssembler< - 'a, - T: RlstScalar, - I: BoundaryIntegrand, - G: GeometryMap, - K: KernelEvaluator, -> { - integrand: &'a I, - kernel: &'a K, - test_evaluator: G, - trial_evaluator: G, - test_table: &'a RlstArray, - trial_table: &'a RlstArray, - k: RlstArray, - test_mapped_pts: RlstArray, - trial_mapped_pts: RlstArray, - test_normals: RlstArray, - trial_normals: RlstArray, - test_jacobians: RlstArray, - trial_jacobians: RlstArray, - test_jdet: Vec, - trial_jdet: Vec, -} - -impl< - 'a, - T: RlstScalar, - I: BoundaryIntegrand, - G: GeometryMap, - K: KernelEvaluator, - > SingularCellPairAssembler<'a, T, I, G, K> -{ - pub fn new( - npts: usize, - deriv_size: usize, - integrand: &'a I, - kernel: &'a K, - test_evaluator: G, - trial_evaluator: G, - test_table: &'a RlstArray, - trial_table: &'a RlstArray, - ) -> Self { - Self { - integrand, - kernel, - test_evaluator, - trial_evaluator, - test_table, - trial_table, - k: rlst_dynamic_array2!(T, [deriv_size, npts]), - test_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), - trial_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), - test_normals: rlst_dynamic_array2!(T::Real, [3, npts]), - trial_normals: rlst_dynamic_array2!(T::Real, [3, npts]), - test_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), - trial_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), - test_jdet: vec![T::Real::zero(); npts], - trial_jdet: vec![T::Real::zero(); npts], - } - } -} -impl< - 'a, - T: RlstScalar, - I: BoundaryIntegrand, - G: GeometryMap, - K: KernelEvaluator, - > SingularCellPairAssembler<'a, T, I, G, K> -{ - // TODO make this a trait - fn assemble( - &mut self, - local_mat: &mut RlstArray, - weights: &[T::Real], - test_cell: usize, - trial_cell: usize, - ) { - self.test_evaluator - .points(test_cell, self.test_mapped_pts.data_mut()); - self.test_evaluator.jacobians_dets_normals( - test_cell, - self.test_jacobians.data_mut(), - &mut self.test_jdet, - self.test_normals.data_mut(), - ); - - self.trial_evaluator - .points(trial_cell, self.trial_mapped_pts.data_mut()); - self.trial_evaluator.jacobians_dets_normals( - trial_cell, - self.trial_jacobians.data_mut(), - &mut self.trial_jdet, - self.trial_normals.data_mut(), - ); - - self.kernel.assemble_pairwise_st( - self.test_mapped_pts.data(), - self.trial_mapped_pts.data(), - self.k.data_mut(), - ); - - let test_geometry = AssemblerGeometry::new( - &self.test_mapped_pts, - &self.test_normals, - &self.test_jacobians, - &self.test_jdet, - ); - let trial_geometry = AssemblerGeometry::new( - &self.trial_mapped_pts, - &self.trial_normals, - &self.trial_jacobians, - &self.trial_jdet, - ); - - for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { - for (test_i, entry) in col.iter_mut().enumerate() { - *entry = T::zero(); - for (index, wt) in weights.iter().enumerate() { - unsafe { - *entry += self.integrand.evaluate_singular( - self.test_table, - self.trial_table, - index, - test_i, - trial_i, - &self.k, - &test_geometry, - &trial_geometry, - ) * num::cast::( - *wt * *self.test_jdet.get_unchecked(index) - * *self.trial_jdet.get_unchecked(index), - ) - .unwrap(); - } - } - } - } - } -} /// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells #[allow(clippy::too_many_arguments)] fn assemble_batch_singular< @@ -294,6 +156,7 @@ fn assemble_batch_singular< trial_evaluator, test_table, trial_table, + weights, ); let mut local_mat = rlst_dynamic_array2!( @@ -304,7 +167,9 @@ fn assemble_batch_singular< ] ); for (test_cell, trial_cell) in cell_pairs { - a.assemble(&mut local_mat, weights, *test_cell, *trial_cell); + a.set_test_cell(*test_cell); + a.set_trial_cell(*trial_cell); + a.assemble(&mut local_mat); let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); @@ -358,106 +223,48 @@ fn assemble_batch_nonadjacent< assert_eq!(trial_grid.geometry_dim(), 3); assert_eq!(trial_grid.topology_dim(), 2); - let mut k = rlst_dynamic_array3!(T, [deriv_size, npts_test, npts_trial]); - let zero = num::cast::(0.0).unwrap(); - let mut test_jdet = vec![zero; npts_test]; - let mut test_mapped_pts = rlst_dynamic_array2!(T::Real, [3, npts_test]); - let mut test_normals = rlst_dynamic_array2!(T::Real, [3, npts_test]); - let mut test_jacobians = rlst_dynamic_array2!(T::Real, [6, npts_test]); - let test_evaluator = test_grid.geometry_map(test_cell_type, test_points.data()); let trial_evaluator = trial_grid.geometry_map(trial_cell_type, trial_points.data()); - let mut trial_jdet = vec![vec![zero; npts_trial]; trial_cells.len()]; - let mut trial_mapped_pts = vec![]; - let mut trial_normals = vec![]; - let mut trial_jacobians = vec![]; - for _i in 0..trial_cells.len() { - trial_mapped_pts.push(rlst_dynamic_array2!(T::Real, [3, npts_trial])); - trial_normals.push(rlst_dynamic_array2!(T::Real, [3, npts_trial])); - trial_jacobians.push(rlst_dynamic_array2!(T::Real, [6, npts_trial])); - } - - for (trial_cell_i, trial_cell) in trial_cells.iter().enumerate() { - trial_evaluator.points(*trial_cell, trial_mapped_pts[trial_cell_i].data_mut()); - trial_evaluator.jacobians_dets_normals( - *trial_cell, - trial_jacobians[trial_cell_i].data_mut(), - &mut trial_jdet[trial_cell_i], - trial_normals[trial_cell_i].data_mut(), - ); - } - - let mut sum: T; - let mut trial_integrands = vec![num::cast::(0.0).unwrap(); npts_trial]; + let mut a = NonsingularCellPairAssembler::new( + npts_test, + npts_trial, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + &test_weights, + &trial_weights, + ); - for test_cell in test_cells { - test_evaluator.points(*test_cell, test_mapped_pts.data_mut()); - test_evaluator.jacobians_dets_normals( - *test_cell, - test_jacobians.data_mut(), - &mut test_jdet, - test_normals.data_mut(), - ); + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); - for (trial_cell_i, trial_cell) in trial_cells.iter().enumerate() { + for (trial_cell_i, trial_cell) in trial_cells.iter().enumerate() { + a.set_test_cell(*trial_cell); + let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); + for test_cell in test_cells { if neighbours(test_grid, trial_grid, *test_cell, *trial_cell) { continue; } - assembler.kernel().assemble_st( - test_mapped_pts.data(), - trial_mapped_pts[trial_cell_i].data(), - k.data_mut(), - ); + a.set_test_cell(*test_cell); + a.assemble(&mut local_mat); let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); - let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - - let test_geometry = AssemblerGeometry::new( - &test_mapped_pts, - &test_normals, - &test_jacobians, - &test_jdet, - ); - let trial_geometry = AssemblerGeometry::new( - &trial_mapped_pts[trial_cell_i], - &trial_normals[trial_cell_i], - &trial_jacobians[trial_cell_i], - &trial_jdet[trial_cell_i], - ); - - for (test_i, test_dof) in test_dofs.iter().enumerate() { - for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { - for (trial_index, trial_wt) in trial_weights.iter().enumerate() { - trial_integrands[trial_index] = num::cast::( - *trial_wt * trial_jdet[trial_cell_i][trial_index], - ) - .unwrap(); - } - sum = num::cast::(0.0).unwrap(); - for (test_index, test_wt) in test_weights.iter().enumerate() { - let test_integrand = - num::cast::(*test_wt * test_jdet[test_index]).unwrap(); - for trial_index in 0..npts_trial { - sum += unsafe { - assembler.integrand().evaluate_nonsingular( - test_table, - trial_table, - test_index, - trial_index, - test_i, - trial_i, - &k, - &test_geometry, - &trial_geometry, - ) * test_integrand - * *trial_integrands.get_unchecked(trial_index) - }; - } - } + + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { unsafe { - *output.data.add(*test_dof + output.shape[0] * *trial_dof) += sum; + *output.data.add(*test_dof + output.shape[0] * *trial_dof) += entry; } } } diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs new file mode 100644 index 00000000..d4323069 --- /dev/null +++ b/src/assembly/boundary/cell_pair_assemblers.rs @@ -0,0 +1,315 @@ +//! Assemblers that assemble the contributions to the global matrix due to a single pair of cells + +use crate::assembly::common::{AssemblerGeometry, RlstArray}; +use crate::traits::{BoundaryIntegrand, CellPairAssembler, KernelEvaluator}; +use ndgrid::traits::GeometryMap; +use num::Zero; +use rlst::{ + rlst_dynamic_array2, rlst_dynamic_array3, DefaultIteratorMut, RawAccess, RawAccessMut, + RlstScalar, +}; + +/// Assembler for the contributions from pairs of neighbouring cells +pub struct SingularCellPairAssembler< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, +> { + integrand: &'a I, + kernel: &'a K, + test_evaluator: G, + trial_evaluator: G, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + k: RlstArray, + test_mapped_pts: RlstArray, + trial_mapped_pts: RlstArray, + test_normals: RlstArray, + trial_normals: RlstArray, + test_jacobians: RlstArray, + trial_jacobians: RlstArray, + test_jdet: Vec, + trial_jdet: Vec, + weights: &'a [T::Real], + test_cell: usize, + trial_cell: usize, +} + +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, + > SingularCellPairAssembler<'a, T, I, G, K> +{ + /// Create new + pub fn new( + npts: usize, + deriv_size: usize, + integrand: &'a I, + kernel: &'a K, + test_evaluator: G, + trial_evaluator: G, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + weights: &'a [T::Real], + ) -> Self { + Self { + integrand, + kernel, + test_evaluator, + trial_evaluator, + test_table, + trial_table, + k: rlst_dynamic_array2!(T, [deriv_size, npts]), + test_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), + trial_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts]), + test_normals: rlst_dynamic_array2!(T::Real, [3, npts]), + trial_normals: rlst_dynamic_array2!(T::Real, [3, npts]), + test_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), + trial_jacobians: rlst_dynamic_array2!(T::Real, [6, npts]), + test_jdet: vec![T::Real::zero(); npts], + trial_jdet: vec![T::Real::zero(); npts], + weights, + test_cell: 0, + trial_cell: 0, + } + } +} +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + G: GeometryMap, + K: KernelEvaluator, + > CellPairAssembler for SingularCellPairAssembler<'a, T, I, G, K> +{ + fn set_test_cell(&mut self, test_cell: usize) { + self.test_cell = test_cell; + self.test_evaluator + .points(test_cell, self.test_mapped_pts.data_mut()); + self.test_evaluator.jacobians_dets_normals( + test_cell, + self.test_jacobians.data_mut(), + &mut self.test_jdet, + self.test_normals.data_mut(), + ); + } + fn set_trial_cell(&mut self, trial_cell: usize) { + self.trial_cell = trial_cell; + self.trial_evaluator + .points(trial_cell, self.trial_mapped_pts.data_mut()); + self.trial_evaluator.jacobians_dets_normals( + trial_cell, + self.trial_jacobians.data_mut(), + &mut self.trial_jdet, + self.trial_normals.data_mut(), + ); + } + fn assemble(&mut self, local_mat: &mut RlstArray) { + self.kernel.assemble_pairwise_st( + self.test_mapped_pts.data(), + self.trial_mapped_pts.data(), + self.k.data_mut(), + ); + + let test_geometry = AssemblerGeometry::new( + &self.test_mapped_pts, + &self.test_normals, + &self.test_jacobians, + &self.test_jdet, + ); + let trial_geometry = AssemblerGeometry::new( + &self.trial_mapped_pts, + &self.trial_normals, + &self.trial_jacobians, + &self.trial_jdet, + ); + + for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { + for (test_i, entry) in col.iter_mut().enumerate() { + *entry = T::zero(); + for (index, wt) in self.weights.iter().enumerate() { + unsafe { + *entry += self.integrand.evaluate_singular( + self.test_table, + self.trial_table, + index, + test_i, + trial_i, + &self.k, + &test_geometry, + &trial_geometry, + ) * num::cast::( + *wt * *self.test_jdet.get_unchecked(index) + * *self.trial_jdet.get_unchecked(index), + ) + .unwrap(); + } + } + } + } + } +} + +/// Assembler for the contributions from pairs of non-neighbouring cells +pub struct NonsingularCellPairAssembler< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, +> { + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + k: RlstArray, + test_mapped_pts: RlstArray, + trial_mapped_pts: RlstArray, + test_normals: RlstArray, + trial_normals: RlstArray, + test_jacobians: RlstArray, + trial_jacobians: RlstArray, + test_jdet: Vec, + trial_jdet: Vec, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + test_cell: usize, + trial_cell: usize, +} + +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > NonsingularCellPairAssembler<'a, T, I, TestG, TrialG, K> +{ + /// Create new + pub fn new( + npts_test: usize, + npts_trial: usize, + deriv_size: usize, + integrand: &'a I, + kernel: &'a K, + test_evaluator: TestG, + trial_evaluator: TrialG, + test_table: &'a RlstArray, + trial_table: &'a RlstArray, + test_weights: &'a [T::Real], + trial_weights: &'a [T::Real], + ) -> Self { + Self { + integrand, + kernel, + test_evaluator, + trial_evaluator, + test_table, + trial_table, + k: rlst_dynamic_array3!(T, [deriv_size, npts_test, npts_trial]), + test_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts_test]), + trial_mapped_pts: rlst_dynamic_array2!(T::Real, [3, npts_trial]), + test_normals: rlst_dynamic_array2!(T::Real, [3, npts_test]), + trial_normals: rlst_dynamic_array2!(T::Real, [3, npts_trial]), + test_jacobians: rlst_dynamic_array2!(T::Real, [6, npts_test]), + trial_jacobians: rlst_dynamic_array2!(T::Real, [6, npts_trial]), + test_jdet: vec![T::Real::zero(); npts_test], + trial_jdet: vec![T::Real::zero(); npts_trial], + test_weights, + trial_weights, + test_cell: 0, + trial_cell: 0, + } + } +} +// TODO: make version of this with trial caching +impl< + 'a, + T: RlstScalar, + I: BoundaryIntegrand, + TestG: GeometryMap, + TrialG: GeometryMap, + K: KernelEvaluator, + > CellPairAssembler for NonsingularCellPairAssembler<'a, T, I, TestG, TrialG, K> +{ + fn set_test_cell(&mut self, test_cell: usize) { + self.test_cell = test_cell; + self.test_evaluator + .points(test_cell, self.test_mapped_pts.data_mut()); + self.test_evaluator.jacobians_dets_normals( + test_cell, + self.test_jacobians.data_mut(), + &mut self.test_jdet, + self.test_normals.data_mut(), + ); + } + fn set_trial_cell(&mut self, trial_cell: usize) { + self.trial_cell = trial_cell; + self.trial_evaluator + .points(trial_cell, self.trial_mapped_pts.data_mut()); + self.trial_evaluator.jacobians_dets_normals( + trial_cell, + self.trial_jacobians.data_mut(), + &mut self.trial_jdet, + self.trial_normals.data_mut(), + ); + } + fn assemble(&mut self, local_mat: &mut RlstArray) { + self.kernel.assemble_st( + self.test_mapped_pts.data(), + self.trial_mapped_pts.data(), + self.k.data_mut(), + ); + + let test_geometry = AssemblerGeometry::new( + &self.test_mapped_pts, + &self.test_normals, + &self.test_jacobians, + &self.test_jdet, + ); + let trial_geometry = AssemblerGeometry::new( + &self.trial_mapped_pts, + &self.trial_normals, + &self.trial_jacobians, + &self.trial_jdet, + ); + + for (trial_i, mut col) in local_mat.col_iter_mut().enumerate() { + for (test_i, entry) in col.iter_mut().enumerate() { + *entry = T::zero(); + for (test_index, test_wt) in self.test_weights.iter().enumerate() { + for (trial_index, trial_wt) in self.trial_weights.iter().enumerate() { + *entry += unsafe { + self.integrand.evaluate_nonsingular( + self.test_table, + self.trial_table, + test_index, + trial_index, + test_i, + trial_i, + &self.k, + &test_geometry, + &trial_geometry, + ) * num::cast::( + *test_wt + * self.test_jdet[test_index] + * *trial_wt + * self.trial_jdet[trial_index], + ) + .unwrap() + }; + } + } + } + } + } +} diff --git a/src/assembly/kernels.rs b/src/assembly/kernels.rs index a4337d5f..c263e51a 100644 --- a/src/assembly/kernels.rs +++ b/src/assembly/kernels.rs @@ -4,22 +4,26 @@ use crate::traits::KernelEvaluator as KernelEvaluatorTrait; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use rlst::RlstScalar; +/// Kernel evaluator pub struct KernelEvaluator> { pub(crate) kernel: K, eval_type: GreenKernelEvalType, } impl> KernelEvaluator { + /// Create new pub fn new(kernel: K, eval_type: GreenKernelEvalType) -> Self { Self { kernel, eval_type } } } impl KernelEvaluator> { + /// Create new Laplace kernel pub fn new_laplace(eval_type: GreenKernelEvalType) -> Self { Self::new(Laplace3dKernel::::new(), eval_type) } } impl> KernelEvaluator> { + /// Create new Helmholtz kernel pub fn new_helmholtz(wavenumber: T::Real, eval_type: GreenKernelEvalType) -> Self { Self::new(Helmholtz3dKernel::::new(wavenumber), eval_type) } diff --git a/src/traits.rs b/src/traits.rs index 30a4ec63..5e9c1b3e 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -3,7 +3,7 @@ mod assembly; mod function; -pub use assembly::{BoundaryIntegrand, CellGeometry, KernelEvaluator}; +pub use assembly::{BoundaryIntegrand, CellGeometry, CellPairAssembler, KernelEvaluator}; pub use function::FunctionSpace; #[cfg(feature = "mpi")] pub use function::ParallelFunctionSpace; diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs index 5fa6440e..2a4cdbb9 100644 --- a/src/traits/assembly.rs +++ b/src/traits/assembly.rs @@ -74,3 +74,13 @@ pub trait KernelEvaluator { result: &mut [Self::T], ); } + +pub trait CellPairAssembler { + //! Assembler for the contributions from a pair of cells + /// Assemble contributions into `local_mat` + fn assemble(&mut self, local_mat: &mut RlstArray); + /// Set the test cell + fn set_test_cell(&mut self, test_cell: usize); + /// Set the trial cell + fn set_trial_cell(&mut self, trial_cell: usize); +} From 0c6be909c5192bcbef70cd26d28be11700b1daf3 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 09:43:08 +0100 Subject: [PATCH 04/13] singular correction --- src/assembly/boundary.rs | 113 ++++++++++++--------------------------- 1 file changed, 33 insertions(+), 80 deletions(-) diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index dcae5613..d35f9db1 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -11,9 +11,7 @@ pub use double_layer::DoubleLayerAssembler; pub use hypersingular::HypersingularAssembler; pub use single_layer::SingleLayerAssembler; -use crate::assembly::common::{ - equal_grids, AssemblerGeometry, RawData2D, RlstArray, SparseMatrixData, -}; +use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData}; use crate::quadrature::duffy::{ quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy, }; @@ -29,12 +27,12 @@ use mpi::traits::Communicator; use ndelement::reference_cell; use ndelement::traits::FiniteElement; use ndelement::types::ReferenceCellType; -use ndgrid::traits::{Entity, GeometryMap, Grid, Topology}; +use ndgrid::traits::{Entity, Grid, Topology}; use ndgrid::types::Ownership; use rayon::prelude::*; use rlst::{ - rlst_dynamic_array2, rlst_dynamic_array3, rlst_dynamic_array4, CsrMatrix, DefaultIterator, - MatrixInverse, RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, + rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, MatrixInverse, + RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, }; use std::collections::HashMap; @@ -248,7 +246,7 @@ fn assemble_batch_nonadjacent< ] ); - for (trial_cell_i, trial_cell) in trial_cells.iter().enumerate() { + for trial_cell in trial_cells { a.set_test_cell(*trial_cell); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); for test_cell in test_cells { @@ -311,90 +309,45 @@ fn assemble_batch_singular_correction< assert_eq!(grid.geometry_dim(), 3); assert_eq!(grid.topology_dim(), 2); - let mut k = rlst_dynamic_array3!(T, [npts_test, deriv_size, npts_trial]); - - let zero = num::cast::(0.0).unwrap(); - - let mut test_jdet = vec![zero; npts_test]; - let mut test_mapped_pts = rlst_dynamic_array2!(T::Real, [3, npts_test]); - let mut test_normals = rlst_dynamic_array2!(T::Real, [3, npts_test]); - let mut test_jacobians = rlst_dynamic_array2!(T::Real, [6, npts_test]); - - let mut trial_jdet = vec![zero; npts_trial]; - let mut trial_mapped_pts = rlst_dynamic_array2!(T::Real, [3, npts_trial]); - let mut trial_normals = rlst_dynamic_array2!(T::Real, [3, npts_trial]); - let mut trial_jacobians = rlst_dynamic_array2!(T::Real, [6, npts_trial]); - let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); - let mut sum: T; - let mut trial_integrands = vec![num::cast::(0.0).unwrap(); npts_trial]; - - for (test_cell, trial_cell) in cell_pairs { - test_evaluator.points(*test_cell, test_mapped_pts.data_mut()); - test_evaluator.jacobians_dets_normals( - *test_cell, - test_jacobians.data_mut(), - &mut test_jdet, - test_normals.data_mut(), - ); + let mut a = NonsingularCellPairAssembler::new( + npts_test, + npts_trial, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + &test_weights, + &trial_weights, + ); - trial_evaluator.points(*trial_cell, trial_mapped_pts.data_mut()); - trial_evaluator.jacobians_dets_normals( - *trial_cell, - trial_jacobians.data_mut(), - &mut trial_jdet, - trial_normals.data_mut(), - ); + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); - assembler.kernel().assemble_st( - test_mapped_pts.data(), - trial_mapped_pts.data(), - k.data_mut(), - ); + for (test_cell, trial_cell) in cell_pairs { + a.set_test_cell(*test_cell); + a.set_test_cell(*trial_cell); - let test_geometry = - AssemblerGeometry::new(&test_mapped_pts, &test_normals, &test_jacobians, &test_jdet); - let trial_geometry = AssemblerGeometry::new( - &trial_mapped_pts, - &trial_normals, - &trial_jacobians, - &trial_jdet, - ); + a.assemble(&mut local_mat); let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - for (test_i, test_dof) in test_dofs.iter().enumerate() { - for (trial_i, trial_dof) in trial_dofs.iter().enumerate() { - for (trial_index, trial_wt) in trial_weights.iter().enumerate() { - trial_integrands[trial_index] = - num::cast::(*trial_wt * trial_jdet[trial_index]).unwrap(); - } - sum = num::cast::(0.0).unwrap(); - for (test_index, test_wt) in test_weights.iter().enumerate() { - let test_integrand = - num::cast::(*test_wt * test_jdet[test_index]).unwrap(); - for trial_index in 0..npts_trial { - sum += unsafe { - assembler.integrand().evaluate_nonsingular( - test_table, - trial_table, - test_index, - trial_index, - test_i, - trial_i, - &k, - &test_geometry, - &trial_geometry, - ) * test_integrand - * *trial_integrands.get_unchecked(trial_index) - }; - } - } + + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { output.rows.push(test_space.global_dof_index(*test_dof)); output.cols.push(trial_space.global_dof_index(*trial_dof)); - output.data.push(sum); + output.data.push(entry); } } } From 57321c7b83da8982c1682ebae66641b07dc85d1c Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 11:53:20 +0100 Subject: [PATCH 05/13] Restructuring traits --- examples/assembly.rs | 2 +- examples/test_parallel_assembly.rs | 2 +- src/assembly.rs | 2 +- src/assembly/boundary.rs | 1119 +--------------- src/assembly/boundary/assemblers.rs | 1141 +++++++++++++++++ .../{ => assemblers}/adjoint_double_layer.rs | 0 .../boundary/{ => assemblers}/double_layer.rs | 0 .../{ => assemblers}/hypersingular.rs | 0 .../boundary/{ => assemblers}/single_layer.rs | 0 src/assembly/boundary/cell_pair_assemblers.rs | 6 +- src/traits.rs | 6 +- src/traits/assembly.rs | 130 +- tests/compare_to_bempp_cl.rs | 2 +- 13 files changed, 1288 insertions(+), 1122 deletions(-) create mode 100644 src/assembly/boundary/assemblers.rs rename src/assembly/boundary/{ => assemblers}/adjoint_double_layer.rs (100%) rename src/assembly/boundary/{ => assemblers}/double_layer.rs (100%) rename src/assembly/boundary/{ => assemblers}/hypersingular.rs (100%) rename src/assembly/boundary/{ => assemblers}/single_layer.rs (100%) diff --git a/examples/assembly.rs b/examples/assembly.rs index 4dd25635..570a4088 100644 --- a/examples/assembly.rs +++ b/examples/assembly.rs @@ -1,6 +1,6 @@ use bempp::assembly::{boundary, boundary::BoundaryAssembler}; use bempp::function::SerialFunctionSpace; -use bempp::traits::FunctionSpace; +use bempp::traits::{BoundaryAssembly, FunctionSpace}; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::{Continuity, ReferenceCellType}; use ndgrid::shapes::regular_sphere; diff --git a/examples/test_parallel_assembly.rs b/examples/test_parallel_assembly.rs index 7f15660a..c2bea395 100644 --- a/examples/test_parallel_assembly.rs +++ b/examples/test_parallel_assembly.rs @@ -7,7 +7,7 @@ use bempp::{ assembly::boundary, assembly::boundary::BoundaryAssembler, function::{ParallelFunctionSpace, SerialFunctionSpace}, - traits::FunctionSpace, + traits::{BoundaryAssembly, FunctionSpace, ParallelBoundaryAssembly}, }; #[cfg(feature = "mpi")] use itertools::izip; diff --git a/src/assembly.rs b/src/assembly.rs index 8a359bbb..941d1cbf 100644 --- a/src/assembly.rs +++ b/src/assembly.rs @@ -10,7 +10,7 @@ mod test { use super::boundary::BoundaryAssembler; use super::*; use crate::function::SerialFunctionSpace; - use crate::traits::FunctionSpace; + use crate::traits::{BoundaryAssembly, FunctionSpace}; use cauchy::{c32, c64}; use ndelement::ciarlet::CiarletElement; use ndelement::ciarlet::LagrangeElementFamily; diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index d35f9db1..922a56d5 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -1,1125 +1,20 @@ //! Assembly of boundary operators -pub(crate) mod adjoint_double_layer; -mod cell_pair_assemblers; -pub(crate) mod double_layer; -pub(crate) mod hypersingular; +mod assemblers; +pub(crate) mod cell_pair_assemblers; mod integrands; -pub(crate) mod single_layer; -pub use adjoint_double_layer::AdjointDoubleLayerAssembler; -pub use double_layer::DoubleLayerAssembler; -pub use hypersingular::HypersingularAssembler; -pub use single_layer::SingleLayerAssembler; - -use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData}; -use crate::quadrature::duffy::{ - quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy, -}; -use crate::quadrature::simplex_rules::simplex_rule; -use crate::quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; -#[cfg(feature = "mpi")] -use crate::traits::ParallelFunctionSpace; -use crate::traits::{BoundaryIntegrand, CellPairAssembler, FunctionSpace, KernelEvaluator}; -use cell_pair_assemblers::{NonsingularCellPairAssembler, SingularCellPairAssembler}; -use itertools::izip; +pub use assemblers::BoundaryAssembler; #[cfg(feature = "mpi")] -use mpi::traits::Communicator; -use ndelement::reference_cell; -use ndelement::traits::FiniteElement; -use ndelement::types::ReferenceCellType; -use ndgrid::traits::{Entity, Grid, Topology}; -use ndgrid::types::Ownership; -use rayon::prelude::*; -use rlst::{ - rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, MatrixInverse, - RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, +pub use assemblers::ParallelBoundaryAssembler; +pub use assemblers::{ + AdjointDoubleLayerAssembler, DoubleLayerAssembler, HypersingularAssembler, SingleLayerAssembler, }; -use std::collections::HashMap; - -fn neighbours( - test_grid: &TestGrid, - trial_grid: &TrialGrid, - test_cell: usize, - trial_cell: usize, -) -> bool { - if !equal_grids(test_grid, trial_grid) { - false - } else { - let test_vertices = trial_grid - .entity(2, test_cell) - .unwrap() - .topology() - .sub_entity_iter(0) - .collect::>(); - for v in trial_grid - .entity(2, trial_cell) - .unwrap() - .topology() - .sub_entity_iter(0) - { - if test_vertices.contains(&v) { - return true; - } - } - false - } -} - -fn get_singular_quadrature_rule( - test_celltype: ReferenceCellType, - trial_celltype: ReferenceCellType, - pairs: &[(usize, usize)], - npoints: usize, -) -> TestTrialNumericalQuadratureDefinition { - if pairs.is_empty() { - panic!("Non-singular rule requested."); - } - let con = CellToCellConnectivity { - connectivity_dimension: match pairs.len() { - 1 => 0, - 2 => 1, - _ => 2, - }, - local_indices: pairs.to_vec(), - }; - match test_celltype { - ReferenceCellType::Triangle => match trial_celltype { - ReferenceCellType::Triangle => triangle_duffy(&con, npoints).unwrap(), - ReferenceCellType::Quadrilateral => { - triangle_quadrilateral_duffy(&con, npoints).unwrap() - } - _ => { - unimplemented!("Only triangles and quadrilaterals are currently supported"); - } - }, - ReferenceCellType::Quadrilateral => match trial_celltype { - ReferenceCellType::Triangle => quadrilateral_triangle_duffy(&con, npoints).unwrap(), - ReferenceCellType::Quadrilateral => quadrilateral_duffy(&con, npoints).unwrap(), - _ => { - unimplemented!("Only triangles and quadrilaterals are currently supported"); - } - }, - _ => { - unimplemented!("Only triangles and quadrilaterals are currently supported"); - } - } -} - -/// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells -#[allow(clippy::too_many_arguments)] -fn assemble_batch_singular< - T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, ->( - assembler: &impl BoundaryAssembler, - deriv_size: usize, - shape: [usize; 2], - trial_cell_type: ReferenceCellType, - test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, - test_space: &impl FunctionSpace, - cell_pairs: &[(usize, usize)], - trial_points: &RlstArray, - test_points: &RlstArray, - weights: &[T::Real], - trial_table: &RlstArray, - test_table: &RlstArray, -) -> SparseMatrixData { - let mut output = SparseMatrixData::::new_known_size( - shape, - cell_pairs.len() - * trial_space.element(trial_cell_type).dim() - * test_space.element(test_cell_type).dim(), - ); - let npts = weights.len(); - debug_assert!(weights.len() == npts); - debug_assert!(test_points.shape()[1] == npts); - debug_assert!(trial_points.shape()[1] == npts); - - let grid = test_space.grid(); - assert_eq!(grid.geometry_dim(), 3); - assert_eq!(grid.topology_dim(), 2); - - let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); - let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); - - let mut a = SingularCellPairAssembler::new( - npts, - deriv_size, - assembler.integrand(), - assembler.kernel(), - test_evaluator, - trial_evaluator, - test_table, - trial_table, - weights, - ); - - let mut local_mat = rlst_dynamic_array2!( - T, - [ - test_space.element(test_cell_type).dim(), - trial_space.element(trial_cell_type).dim() - ] - ); - for (test_cell, trial_cell) in cell_pairs { - a.set_test_cell(*test_cell); - a.set_trial_cell(*trial_cell); - a.assemble(&mut local_mat); - - let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); - let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - - for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { - for (test_dof, entry) in izip!(test_dofs, col.iter()) { - output.rows.push(test_space.global_dof_index(*test_dof)); - output.cols.push(trial_space.global_dof_index(*trial_dof)); - output.data.push(entry); - } - } - } - - output -} - -/// Assemble the contribution to the terms of a matrix for a batch of non-adjacent cells -#[allow(clippy::too_many_arguments)] -fn assemble_batch_nonadjacent< - T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, ->( - assembler: &impl BoundaryAssembler, - deriv_size: usize, - output: &RawData2D, - trial_cell_type: ReferenceCellType, - test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, - trial_cells: &[usize], - test_space: &impl FunctionSpace, - test_cells: &[usize], - trial_points: &RlstArray, - trial_weights: &[T::Real], - test_points: &RlstArray, - test_weights: &[T::Real], - trial_table: &RlstArray, - test_table: &RlstArray, -) -> usize { - let npts_test = test_weights.len(); - let npts_trial = trial_weights.len(); - debug_assert!(test_points.shape()[1] == npts_test); - debug_assert!(trial_points.shape()[1] == npts_trial); - - let test_grid = test_space.grid(); - let trial_grid = trial_space.grid(); - - assert_eq!(test_grid.geometry_dim(), 3); - assert_eq!(test_grid.topology_dim(), 2); - assert_eq!(trial_grid.geometry_dim(), 3); - assert_eq!(trial_grid.topology_dim(), 2); - - let test_evaluator = test_grid.geometry_map(test_cell_type, test_points.data()); - let trial_evaluator = trial_grid.geometry_map(trial_cell_type, trial_points.data()); - - let mut a = NonsingularCellPairAssembler::new( - npts_test, - npts_trial, - deriv_size, - assembler.integrand(), - assembler.kernel(), - test_evaluator, - trial_evaluator, - test_table, - trial_table, - &test_weights, - &trial_weights, - ); - - let mut local_mat = rlst_dynamic_array2!( - T, - [ - test_space.element(test_cell_type).dim(), - trial_space.element(trial_cell_type).dim() - ] - ); - - for trial_cell in trial_cells { - a.set_test_cell(*trial_cell); - let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - for test_cell in test_cells { - if neighbours(test_grid, trial_grid, *test_cell, *trial_cell) { - continue; - } - - a.set_test_cell(*test_cell); - a.assemble(&mut local_mat); - - let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); - - for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { - for (test_dof, entry) in izip!(test_dofs, col.iter()) { - unsafe { - *output.data.add(*test_dof + output.shape[0] * *trial_dof) += entry; - } - } - } - } - } - 1 -} - -/// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells if an (incorrect) non-singular quadrature rule was used -#[allow(clippy::too_many_arguments)] -fn assemble_batch_singular_correction< - T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, ->( - assembler: &impl BoundaryAssembler, - deriv_size: usize, - shape: [usize; 2], - trial_cell_type: ReferenceCellType, - test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, - test_space: &impl FunctionSpace, - cell_pairs: &[(usize, usize)], - trial_points: &RlstArray, - trial_weights: &[T::Real], - test_points: &RlstArray, - test_weights: &[T::Real], - trial_table: &RlstArray, - test_table: &RlstArray, -) -> SparseMatrixData { - let mut output = SparseMatrixData::::new_known_size( - shape, - cell_pairs.len() - * trial_space.element(trial_cell_type).dim() - * test_space.element(test_cell_type).dim(), - ); - let npts_test = test_weights.len(); - let npts_trial = trial_weights.len(); - debug_assert!(test_points.shape()[1] == npts_test); - debug_assert!(trial_points.shape()[1] == npts_trial); - - let grid = test_space.grid(); - assert_eq!(grid.geometry_dim(), 3); - assert_eq!(grid.topology_dim(), 2); - - let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); - let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); - - let mut a = NonsingularCellPairAssembler::new( - npts_test, - npts_trial, - deriv_size, - assembler.integrand(), - assembler.kernel(), - test_evaluator, - trial_evaluator, - test_table, - trial_table, - &test_weights, - &trial_weights, - ); - - let mut local_mat = rlst_dynamic_array2!( - T, - [ - test_space.element(test_cell_type).dim(), - trial_space.element(trial_cell_type).dim() - ] - ); - - for (test_cell, trial_cell) in cell_pairs { - a.set_test_cell(*test_cell); - a.set_test_cell(*trial_cell); - - a.assemble(&mut local_mat); - - let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); - let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); - - for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { - for (test_dof, entry) in izip!(test_dofs, col.iter()) { - output.rows.push(test_space.global_dof_index(*test_dof)); - output.cols.push(trial_space.global_dof_index(*trial_dof)); - output.data.push(entry); - } - } - } - output -} - -fn get_pairs_if_smallest( - test_cell: &impl Entity, - trial_cell: &impl Entity, - vertex: usize, -) -> Option> { - let mut pairs = vec![]; - for (trial_i, trial_v) in trial_cell.topology().sub_entity_iter(0).enumerate() { - for (test_i, test_v) in test_cell.topology().sub_entity_iter(0).enumerate() { - if test_v == trial_v { - if test_v < vertex { - return None; - } - pairs.push((test_i, trial_i)); - } - } - } - Some(pairs) -} - -/// Options for a boundary assembler -pub struct BoundaryAssemblerOptions { - /// Number of points used in quadrature for non-singular integrals - quadrature_degrees: HashMap, - /// Quadrature degrees to be used for singular integrals - singular_quadrature_degrees: HashMap<(ReferenceCellType, ReferenceCellType), usize>, - /// Maximum size of each batch of cells to send to an assembly function - batch_size: usize, -} - -impl Default for BoundaryAssemblerOptions { - fn default() -> Self { - use ReferenceCellType::{Quadrilateral, Triangle}; - Self { - quadrature_degrees: HashMap::from([(Triangle, 37), (Quadrilateral, 37)]), - singular_quadrature_degrees: HashMap::from([ - ((Triangle, Triangle), 4), - ((Quadrilateral, Quadrilateral), 4), - ((Quadrilateral, Triangle), 4), - ((Triangle, Quadrilateral), 4), - ]), - batch_size: 128, - } - } -} - -pub trait BoundaryAssembler: Sync + Sized { - //! Boundary assembler - //! - //! Assemble operators by processing batches of cells in parallel - - /// Scalar type - type T: RlstScalar + MatrixInverse; - /// Integrand type - type Integrand: BoundaryIntegrand; - /// Kernel type - type Kernel: KernelEvaluator; - /// Number of derivatives - const DERIV_SIZE: usize; - /// Number of derivatives needed in basis function tables - const TABLE_DERIVS: usize; - - /// Get integrand - fn integrand(&self) -> &Self::Integrand; - - /// Get integrand - fn kernel(&self) -> &Self::Kernel; - - /// Get assembler options - fn options(&self) -> &BoundaryAssemblerOptions; - - /// Get mutable assembler options - fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions; - - /// Set (non-singular) quadrature degree for a cell type - fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) { - *self - .options_mut() - .quadrature_degrees - .get_mut(&cell) - .unwrap() = degree; - } - - /// Set singular quadrature degree for a pair of cell types - fn singular_quadrature_degree( - &mut self, - cells: (ReferenceCellType, ReferenceCellType), - degree: usize, - ) { - *self - .options_mut() - .singular_quadrature_degrees - .get_mut(&cells) - .unwrap() = degree; - } - - /// Set the maximum size of a batch of cells to send to an assembly function - fn batch_size(&mut self, size: usize) { - self.options_mut().batch_size = size; - } - - /// Assemble the singular contributions - fn assemble_singular< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> SparseMatrixData { - if !equal_grids(test_space.grid(), trial_space.grid()) { - // If the test and trial grids are different, there are no neighbouring triangles - return SparseMatrixData::new(shape); - } - - if !trial_space.is_serial() || !test_space.is_serial() { - panic!("Dense assembly can only be used for function spaces stored in serial"); - } - if shape[0] != test_space.global_size() || shape[1] != trial_space.global_size() { - panic!("Matrix has wrong shape"); - } - - let grid = test_space.grid(); - - let mut qweights = vec![]; - let mut trial_points = vec![]; - let mut test_points = vec![]; - let mut trial_tables = vec![]; - let mut test_tables = vec![]; - let mut test_cell_types = vec![]; - let mut trial_cell_types = vec![]; - - let mut cell_blocks = vec![]; - - let mut pair_indices = HashMap::new(); - - for test_cell_type in grid.entity_types(2) { - for trial_cell_type in grid.entity_types(2) { - let qdegree = self.options().singular_quadrature_degrees - [&(*test_cell_type, *trial_cell_type)]; - let offset = qweights.len(); - - let mut possible_pairs = vec![]; - // Vertex-adjacent - for i in 0..reference_cell::entity_counts(*test_cell_type)[0] { - for j in 0..reference_cell::entity_counts(*trial_cell_type)[0] { - possible_pairs.push(vec![(i, j)]); - } - } - // edge-adjacent - for test_e in reference_cell::edges(*test_cell_type) { - for trial_e in reference_cell::edges(*trial_cell_type) { - possible_pairs.push(vec![(test_e[0], trial_e[0]), (test_e[1], trial_e[1])]); - possible_pairs.push(vec![(test_e[1], trial_e[0]), (test_e[0], trial_e[1])]); - } - } - // Same cell - if test_cell_type == trial_cell_type { - possible_pairs.push( - (0..reference_cell::entity_counts(*test_cell_type)[0]) - .map(&|i| (i, i)) - .collect::>(), - ); - } - - for (i, pairs) in possible_pairs.iter().enumerate() { - pair_indices.insert( - (*test_cell_type, *trial_cell_type, pairs.clone()), - offset + i, - ); - test_cell_types.push(*test_cell_type); - trial_cell_types.push(*trial_cell_type); - } - - for pairs in &possible_pairs { - let qrule = get_singular_quadrature_rule( - *test_cell_type, - *trial_cell_type, - pairs, - qdegree, - ); - let npts = qrule.weights.len(); - - let mut points = rlst_dynamic_array2!(::Real, [2, npts]); - for i in 0..npts { - for j in 0..2 { - *points.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule.trial_points[2 * i + j], - ) - .unwrap(); - } - } - let trial_element = trial_space.element(*trial_cell_type); - let mut table = rlst_dynamic_array4!( - Self::T, - trial_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1]) - ); - trial_element.tabulate(&points, Self::TABLE_DERIVS, &mut table); - trial_points.push(points); - trial_tables.push(table); - - let mut points = rlst_dynamic_array2!(::Real, [2, npts]); - for i in 0..npts { - for j in 0..2 { - *points.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule.test_points[2 * i + j], - ) - .unwrap(); - } - } - let test_element = test_space.element(*test_cell_type); - let mut table = rlst_dynamic_array4!( - Self::T, - test_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1]) - ); - test_element.tabulate(&points, Self::TABLE_DERIVS, &mut table); - test_points.push(points); - test_tables.push(table); - qweights.push( - qrule - .weights - .iter() - .map(|w| num::cast::::Real>(*w).unwrap()) - .collect::>(), - ); - } - } - } - let mut cell_pairs: Vec> = vec![vec![]; pair_indices.len()]; - for vertex in grid.entity_iter(0) { - for test_cell_index in vertex.topology().connected_entity_iter(2) { - let test_cell = grid.entity(2, test_cell_index).unwrap(); - if test_cell.ownership() == Ownership::Owned { - let test_cell_type = test_cell.entity_type(); - for trial_cell_index in vertex.topology().connected_entity_iter(2) { - let trial_cell = grid.entity(2, trial_cell_index).unwrap(); - let trial_cell_type = trial_cell.entity_type(); - if let Some(pairs) = - get_pairs_if_smallest(&test_cell, &trial_cell, vertex.local_index()) - { - cell_pairs[pair_indices[&(test_cell_type, trial_cell_type, pairs)]] - .push((test_cell_index, trial_cell_index)); - } - } - } - } - } - - let batch_size = self.options().batch_size; - for (i, cells) in cell_pairs.iter().enumerate() { - let mut start = 0; - while start < cells.len() { - let end = std::cmp::min(start + batch_size, cells.len()); - cell_blocks.push((i, &cells[start..end])); - start = end; - } - } - cell_blocks - .into_par_iter() - .map(|(i, cell_block)| { - assemble_batch_singular::( - self, - Self::DERIV_SIZE, - shape, - trial_cell_types[i], - test_cell_types[i], - trial_space, - test_space, - cell_block, - &trial_points[i], - &test_points[i], - &qweights[i], - &trial_tables[i], - &test_tables[i], - ) - }) - .reduce( - || SparseMatrixData::::new(shape), - |mut a, b| { - a.add(b); - a - }, - ) - } - - /// Assemble the singular correction - /// - /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> SparseMatrixData { - if !equal_grids(test_space.grid(), trial_space.grid()) { - // If the test and trial grids are different, there are no neighbouring triangles - return SparseMatrixData::new(shape); - } - - if !trial_space.is_serial() || !test_space.is_serial() { - panic!("Dense assembly can only be used for function spaces stored in serial"); - } - if shape[0] != test_space.global_size() || shape[1] != trial_space.global_size() { - panic!("Matrix has wrong shape"); - } - - let grid = test_space.grid(); - - let mut qweights_test = vec![]; - let mut qweights_trial = vec![]; - let mut qpoints_test = vec![]; - let mut qpoints_trial = vec![]; - let mut test_tables = vec![]; - let mut trial_tables = vec![]; - let mut test_cell_types = vec![]; - let mut trial_cell_types = vec![]; - - let mut cell_blocks = vec![]; - - let mut cell_type_indices = HashMap::new(); - - for test_cell_type in grid.entity_types(2) { - let npts_test = self.options().quadrature_degrees[test_cell_type]; - for trial_cell_type in grid.entity_types(2) { - let npts_trial = self.options().quadrature_degrees[trial_cell_type]; - test_cell_types.push(*test_cell_type); - trial_cell_types.push(*trial_cell_type); - cell_type_indices.insert((*test_cell_type, *trial_cell_type), qweights_test.len()); - - let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap(); - let mut test_pts = - rlst_dynamic_array2!(::Real, [2, npts_test]); - for i in 0..npts_test { - for j in 0..2 { - *test_pts.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule_test.points[2 * i + j], - ) - .unwrap(); - } - } - qweights_test.push( - qrule_test - .weights - .iter() - .map(|w| num::cast::::Real>(*w).unwrap()) - .collect::>(), - ); - let test_element = test_space.element(*test_cell_type); - let mut test_table = rlst_dynamic_array4!( - Self::T, - test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test) - ); - test_element.tabulate(&test_pts, Self::TABLE_DERIVS, &mut test_table); - test_tables.push(test_table); - qpoints_test.push(test_pts); - - let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap(); - let mut trial_pts = - rlst_dynamic_array2!(::Real, [2, npts_trial]); - for i in 0..npts_trial { - for j in 0..2 { - *trial_pts.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule_trial.points[2 * i + j], - ) - .unwrap(); - } - } - qweights_trial.push( - qrule_trial - .weights - .iter() - .map(|w| num::cast::::Real>(*w).unwrap()) - .collect::>(), - ); - let trial_element = trial_space.element(*trial_cell_type); - let mut trial_table = rlst_dynamic_array4!( - Self::T, - trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial) - ); - trial_element.tabulate(&trial_pts, Self::TABLE_DERIVS, &mut trial_table); - trial_tables.push(trial_table); - qpoints_trial.push(trial_pts); - } - } - let mut cell_pairs: Vec> = vec![vec![]; qweights_test.len()]; - - for vertex in grid.entity_iter(0) { - for test_cell_index in vertex.topology().connected_entity_iter(2) { - let test_cell = grid.entity(2, test_cell_index).unwrap(); - let test_cell_type = test_cell.entity_type(); - if test_cell.ownership() == Ownership::Owned { - for trial_cell_index in vertex.topology().connected_entity_iter(2) { - let trial_cell = grid.entity(2, trial_cell_index).unwrap(); - let trial_cell_type = trial_cell.entity_type(); - - if get_pairs_if_smallest(&test_cell, &trial_cell, vertex.local_index()) - .is_some() - { - cell_pairs[cell_type_indices[&(test_cell_type, trial_cell_type)]] - .push((test_cell_index, trial_cell_index)); - } - } - } - } - } - let batch_size = self.options().batch_size; - for (i, cells) in cell_pairs.iter().enumerate() { - let mut start = 0; - while start < cells.len() { - let end = std::cmp::min(start + batch_size, cells.len()); - cell_blocks.push((i, &cells[start..end])); - start = end; - } - } - - cell_blocks - .into_par_iter() - .map(|(i, cell_block)| { - assemble_batch_singular_correction::( - self, - Self::DERIV_SIZE, - shape, - trial_cell_types[i], - test_cell_types[i], - trial_space, - test_space, - cell_block, - &qpoints_trial[i], - &qweights_trial[i], - &qpoints_test[i], - &qweights_test[i], - &trial_tables[i], - &test_tables[i], - ) - }) - .reduce( - || SparseMatrixData::::new(shape), - |mut a, b| { - a.add(b); - a - }, - ) - } - - /// Assemble the singular contributions into a dense matrix - fn assemble_singular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) { - let sparse_matrix = self.assemble_singular::( - output.shape(), - trial_space, - test_space, - ); - let data = sparse_matrix.data; - let rows = sparse_matrix.rows; - let cols = sparse_matrix.cols; - for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { - *output.get_mut([*i, *j]).unwrap() += *value; - } - } - - /// Assemble the singular contributions into a CSR sparse matrix - fn assemble_singular_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> CsrMatrix { - let shape = [test_space.global_size(), trial_space.global_size()]; - let sparse_matrix = - self.assemble_singular::(shape, trial_space, test_space); - - CsrMatrix::::from_aij( - sparse_matrix.shape, - &sparse_matrix.rows, - &sparse_matrix.cols, - &sparse_matrix.data, - ) - .unwrap() - } - - #[cfg(feature = "mpi")] - /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers - fn parallel_assemble_singular_into_csr< - 'a, - C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, - >( - &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), - ) -> CsrMatrix { - self.assemble_singular_into_csr(trial_space.local_space(), test_space.local_space()) - } - - /// Assemble the singular correction into a dense matrix - /// - /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) { - let sparse_matrix = self.assemble_singular_correction::( - output.shape(), - trial_space, - test_space, - ); - let data = sparse_matrix.data; - let rows = sparse_matrix.rows; - let cols = sparse_matrix.cols; - for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { - *output.get_mut([*i, *j]).unwrap() += *value; - } - } - - /// Assemble the singular correction into a CSR matrix - /// - /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) -> CsrMatrix { - let shape = [test_space.global_size(), trial_space.global_size()]; - let sparse_matrix = self.assemble_singular_correction::( - shape, - trial_space, - test_space, - ); - - CsrMatrix::::from_aij( - sparse_matrix.shape, - &sparse_matrix.rows, - &sparse_matrix.cols, - &sparse_matrix.data, - ) - .unwrap() - } - - #[cfg(feature = "mpi")] - /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers - fn parallel_assemble_singular_correction_into_csr< - 'a, - C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, - >( - &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), - ) -> CsrMatrix { - self.assemble_singular_correction_into_csr( - trial_space.local_space(), - test_space.local_space(), - ) - } - - /// Assemble into a dense matrix - fn assemble_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - ) { - let test_colouring = test_space.cell_colouring(); - let trial_colouring = trial_space.cell_colouring(); - - self.assemble_nonsingular_into_dense::( - output, - trial_space, - test_space, - &trial_colouring, - &test_colouring, - ); - self.assemble_singular_into_dense::( - output, - trial_space, - test_space, - ); - } - - /// Assemble the non-singular contributions into a dense matrix - fn assemble_nonsingular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( - &self, - output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), - trial_colouring: &HashMap>>, - test_colouring: &HashMap>>, - ) { - if !trial_space.is_serial() || !test_space.is_serial() { - panic!("Dense assembly can only be used for function spaces stored in serial"); - } - if output.shape()[0] != test_space.global_size() - || output.shape()[1] != trial_space.global_size() - { - panic!("Matrix has wrong shape"); - } - - let batch_size = self.options().batch_size; - - for test_cell_type in test_space.grid().entity_types(2) { - let npts_test = self.options().quadrature_degrees[test_cell_type]; - for trial_cell_type in trial_space.grid().entity_types(2) { - let npts_trial = self.options().quadrature_degrees[trial_cell_type]; - let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap(); - let mut qpoints_test = - rlst_dynamic_array2!(::Real, [2, npts_test]); - for i in 0..npts_test { - for j in 0..2 { - *qpoints_test.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule_test.points[2 * i + j], - ) - .unwrap(); - } - } - let qweights_test = qrule_test - .weights - .iter() - .map(|w| num::cast::::Real>(*w).unwrap()) - .collect::>(); - let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap(); - let mut qpoints_trial = - rlst_dynamic_array2!(::Real, [2, npts_trial]); - for i in 0..npts_trial { - for j in 0..2 { - *qpoints_trial.get_mut([j, i]).unwrap() = - num::cast::::Real>( - qrule_trial.points[2 * i + j], - ) - .unwrap(); - } - } - let qweights_trial = qrule_trial - .weights - .iter() - .map(|w| num::cast::::Real>(*w).unwrap()) - .collect::>(); - - let test_element = test_space.element(*test_cell_type); - let mut test_table = rlst_dynamic_array4!( - Self::T, - test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test) - ); - test_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut test_table); - - let trial_element = trial_space.element(*trial_cell_type); - let mut trial_table = rlst_dynamic_array4!( - Self::T, - trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial) - ); - trial_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut trial_table); - - let output_raw = RawData2D { - data: output.data_mut().as_mut_ptr(), - shape: output.shape(), - }; - - for test_c in &test_colouring[test_cell_type] { - for trial_c in &trial_colouring[trial_cell_type] { - let mut test_cells: Vec<&[usize]> = vec![]; - let mut trial_cells: Vec<&[usize]> = vec![]; - - let mut test_start = 0; - while test_start < test_c.len() { - let test_end = if test_start + batch_size < test_c.len() { - test_start + batch_size - } else { - test_c.len() - }; - - let mut trial_start = 0; - while trial_start < trial_c.len() { - let trial_end = if trial_start + batch_size < trial_c.len() { - trial_start + batch_size - } else { - trial_c.len() - }; - test_cells.push(&test_c[test_start..test_end]); - trial_cells.push(&trial_c[trial_start..trial_end]); - trial_start = trial_end; - } - test_start = test_end - } - - let numtasks = test_cells.len(); - let r: usize = (0..numtasks) - .into_par_iter() - .map(&|t| { - assemble_batch_nonadjacent::( - self, - Self::DERIV_SIZE, - &output_raw, - *test_cell_type, - *trial_cell_type, - trial_space, - trial_cells[t], - test_space, - test_cells[t], - &qpoints_trial, - &qweights_trial, - &qpoints_test, - &qweights_test, - &trial_table, - &test_table, - ) - }) - .sum(); - assert_eq!(r, numtasks); - } - } - } - } - } -} #[cfg(test)] mod test { use super::*; use crate::function::SerialFunctionSpace; - use crate::traits::FunctionSpace; + use crate::traits::{BoundaryAssembly, FunctionSpace}; use approx::*; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::Continuity; diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs new file mode 100644 index 00000000..d758bdb6 --- /dev/null +++ b/src/assembly/boundary/assemblers.rs @@ -0,0 +1,1141 @@ +//! Boundary assemblers +pub mod adjoint_double_layer; +pub mod double_layer; +pub mod hypersingular; +pub mod single_layer; +pub use adjoint_double_layer::AdjointDoubleLayerAssembler; +pub use double_layer::DoubleLayerAssembler; +pub use hypersingular::HypersingularAssembler; +pub use single_layer::SingleLayerAssembler; + +use super::cell_pair_assemblers::{NonsingularCellPairAssembler, SingularCellPairAssembler}; +use crate::assembly::common::{equal_grids, RawData2D, RlstArray, SparseMatrixData}; +use crate::quadrature::duffy::{ + quadrilateral_duffy, quadrilateral_triangle_duffy, triangle_duffy, triangle_quadrilateral_duffy, +}; +use crate::quadrature::simplex_rules::simplex_rule; +use crate::quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; +#[cfg(feature = "mpi")] +use crate::traits::ParallelFunctionSpace; +use crate::traits::{BoundaryIntegrand, CellPairAssembler, FunctionSpace, KernelEvaluator}; +use itertools::izip; +#[cfg(feature = "mpi")] +use mpi::traits::Communicator; +use ndelement::reference_cell; +use ndelement::traits::FiniteElement; +use ndelement::types::ReferenceCellType; +use ndgrid::traits::{Entity, Grid, Topology}; +use ndgrid::types::Ownership; +use rayon::prelude::*; +use rlst::{ + rlst_dynamic_array2, rlst_dynamic_array4, CsrMatrix, DefaultIterator, MatrixInverse, + RandomAccessMut, RawAccess, RawAccessMut, RlstScalar, Shape, +}; +use std::collections::HashMap; + +fn neighbours( + test_grid: &TestGrid, + trial_grid: &TrialGrid, + test_cell: usize, + trial_cell: usize, +) -> bool { + if !equal_grids(test_grid, trial_grid) { + false + } else { + let test_vertices = trial_grid + .entity(2, test_cell) + .unwrap() + .topology() + .sub_entity_iter(0) + .collect::>(); + for v in trial_grid + .entity(2, trial_cell) + .unwrap() + .topology() + .sub_entity_iter(0) + { + if test_vertices.contains(&v) { + return true; + } + } + false + } +} + +fn get_singular_quadrature_rule( + test_celltype: ReferenceCellType, + trial_celltype: ReferenceCellType, + pairs: &[(usize, usize)], + npoints: usize, +) -> TestTrialNumericalQuadratureDefinition { + if pairs.is_empty() { + panic!("Non-singular rule requested."); + } + let con = CellToCellConnectivity { + connectivity_dimension: match pairs.len() { + 1 => 0, + 2 => 1, + _ => 2, + }, + local_indices: pairs.to_vec(), + }; + match test_celltype { + ReferenceCellType::Triangle => match trial_celltype { + ReferenceCellType::Triangle => triangle_duffy(&con, npoints).unwrap(), + ReferenceCellType::Quadrilateral => { + triangle_quadrilateral_duffy(&con, npoints).unwrap() + } + _ => { + unimplemented!("Only triangles and quadrilaterals are currently supported"); + } + }, + ReferenceCellType::Quadrilateral => match trial_celltype { + ReferenceCellType::Triangle => quadrilateral_triangle_duffy(&con, npoints).unwrap(), + ReferenceCellType::Quadrilateral => quadrilateral_duffy(&con, npoints).unwrap(), + _ => { + unimplemented!("Only triangles and quadrilaterals are currently supported"); + } + }, + _ => { + unimplemented!("Only triangles and quadrilaterals are currently supported"); + } + } +} + +fn make_cell_blocks( + f: F, + size: usize, + grid: &impl Grid, + batch_size: usize, +) -> Vec<(usize, Vec<(usize, usize)>)> +where + F: Fn(ReferenceCellType, ReferenceCellType, Vec<(usize, usize)>) -> usize, +{ + let mut cell_pairs = vec![vec![]; size]; + + for vertex in grid.entity_iter(0) { + for test_cell_index in vertex.topology().connected_entity_iter(2) { + let test_cell = grid.entity(2, test_cell_index).unwrap(); + let test_cell_type = test_cell.entity_type(); + if test_cell.ownership() == Ownership::Owned { + for trial_cell_index in vertex.topology().connected_entity_iter(2) { + let trial_cell = grid.entity(2, trial_cell_index).unwrap(); + let trial_cell_type = trial_cell.entity_type(); + + if let Some(pairs) = + get_pairs_if_smallest(&test_cell, &trial_cell, vertex.local_index()) + { + cell_pairs[f(test_cell_type, trial_cell_type, pairs)] + .push((test_cell_index, trial_cell_index)); + } + } + } + } + } + let mut cell_blocks = vec![]; + + for (i, cells) in cell_pairs.iter().enumerate() { + let mut start = 0; + while start < cells.len() { + let end = std::cmp::min(start + batch_size, cells.len()); + cell_blocks.push((i, cells[start..end].to_vec())); + start = end; + } + } + + cell_blocks +} + +/// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells +#[allow(clippy::too_many_arguments)] +fn assemble_batch_singular< + T: RlstScalar + MatrixInverse, + TestGrid: Grid, + TrialGrid: Grid, + Element: FiniteElement + Sync, +>( + assembler: &impl BoundaryAssembler, + deriv_size: usize, + shape: [usize; 2], + trial_cell_type: ReferenceCellType, + test_cell_type: ReferenceCellType, + trial_space: &impl FunctionSpace, + test_space: &impl FunctionSpace, + cell_pairs: &[(usize, usize)], + trial_points: &RlstArray, + test_points: &RlstArray, + weights: &[T::Real], + trial_table: &RlstArray, + test_table: &RlstArray, +) -> SparseMatrixData { + let mut output = SparseMatrixData::::new_known_size( + shape, + cell_pairs.len() + * trial_space.element(trial_cell_type).dim() + * test_space.element(test_cell_type).dim(), + ); + let npts = weights.len(); + debug_assert!(weights.len() == npts); + debug_assert!(test_points.shape()[1] == npts); + debug_assert!(trial_points.shape()[1] == npts); + + let grid = test_space.grid(); + assert_eq!(grid.geometry_dim(), 3); + assert_eq!(grid.topology_dim(), 2); + + let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); + let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); + + let mut a = SingularCellPairAssembler::new( + npts, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + weights, + ); + + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); + for (test_cell, trial_cell) in cell_pairs { + a.set_test_cell(*test_cell); + a.set_trial_cell(*trial_cell); + a.assemble(&mut local_mat); + + let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); + let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); + + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { + output.rows.push(test_space.global_dof_index(*test_dof)); + output.cols.push(trial_space.global_dof_index(*trial_dof)); + output.data.push(entry); + } + } + } + + output +} + +/// Assemble the contribution to the terms of a matrix for a batch of non-adjacent cells +#[allow(clippy::too_many_arguments)] +fn assemble_batch_nonadjacent< + T: RlstScalar + MatrixInverse, + TestGrid: Grid, + TrialGrid: Grid, + Element: FiniteElement + Sync, +>( + assembler: &impl BoundaryAssembler, + deriv_size: usize, + output: &RawData2D, + trial_cell_type: ReferenceCellType, + test_cell_type: ReferenceCellType, + trial_space: &impl FunctionSpace, + trial_cells: &[usize], + test_space: &impl FunctionSpace, + test_cells: &[usize], + trial_points: &RlstArray, + trial_weights: &[T::Real], + test_points: &RlstArray, + test_weights: &[T::Real], + trial_table: &RlstArray, + test_table: &RlstArray, +) -> usize { + let npts_test = test_weights.len(); + let npts_trial = trial_weights.len(); + debug_assert!(test_points.shape()[1] == npts_test); + debug_assert!(trial_points.shape()[1] == npts_trial); + + let test_grid = test_space.grid(); + let trial_grid = trial_space.grid(); + + assert_eq!(test_grid.geometry_dim(), 3); + assert_eq!(test_grid.topology_dim(), 2); + assert_eq!(trial_grid.geometry_dim(), 3); + assert_eq!(trial_grid.topology_dim(), 2); + + let test_evaluator = test_grid.geometry_map(test_cell_type, test_points.data()); + let trial_evaluator = trial_grid.geometry_map(trial_cell_type, trial_points.data()); + + let mut a = NonsingularCellPairAssembler::new( + npts_test, + npts_trial, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + &test_weights, + &trial_weights, + ); + + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); + + for trial_cell in trial_cells { + a.set_test_cell(*trial_cell); + let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); + for test_cell in test_cells { + if neighbours(test_grid, trial_grid, *test_cell, *trial_cell) { + continue; + } + + a.set_test_cell(*test_cell); + a.assemble(&mut local_mat); + + let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); + + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { + unsafe { + *output.data.add(*test_dof + output.shape[0] * *trial_dof) += entry; + } + } + } + } + } + 1 +} + +/// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells if an (incorrect) non-singular quadrature rule was used +#[allow(clippy::too_many_arguments)] +fn assemble_batch_singular_correction< + T: RlstScalar + MatrixInverse, + TestGrid: Grid, + TrialGrid: Grid, + Element: FiniteElement + Sync, +>( + assembler: &impl BoundaryAssembler, + deriv_size: usize, + shape: [usize; 2], + trial_cell_type: ReferenceCellType, + test_cell_type: ReferenceCellType, + trial_space: &impl FunctionSpace, + test_space: &impl FunctionSpace, + cell_pairs: &[(usize, usize)], + trial_points: &RlstArray, + trial_weights: &[T::Real], + test_points: &RlstArray, + test_weights: &[T::Real], + trial_table: &RlstArray, + test_table: &RlstArray, +) -> SparseMatrixData { + let mut output = SparseMatrixData::::new_known_size( + shape, + cell_pairs.len() + * trial_space.element(trial_cell_type).dim() + * test_space.element(test_cell_type).dim(), + ); + let npts_test = test_weights.len(); + let npts_trial = trial_weights.len(); + debug_assert!(test_points.shape()[1] == npts_test); + debug_assert!(trial_points.shape()[1] == npts_trial); + + let grid = test_space.grid(); + assert_eq!(grid.geometry_dim(), 3); + assert_eq!(grid.topology_dim(), 2); + + let test_evaluator = grid.geometry_map(test_cell_type, test_points.data()); + let trial_evaluator = grid.geometry_map(trial_cell_type, trial_points.data()); + + let mut a = NonsingularCellPairAssembler::new( + npts_test, + npts_trial, + deriv_size, + assembler.integrand(), + assembler.kernel(), + test_evaluator, + trial_evaluator, + test_table, + trial_table, + &test_weights, + &trial_weights, + ); + + let mut local_mat = rlst_dynamic_array2!( + T, + [ + test_space.element(test_cell_type).dim(), + trial_space.element(trial_cell_type).dim() + ] + ); + + for (test_cell, trial_cell) in cell_pairs { + a.set_test_cell(*test_cell); + a.set_test_cell(*trial_cell); + + a.assemble(&mut local_mat); + + let test_dofs = test_space.cell_dofs(*test_cell).unwrap(); + let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); + + for (trial_dof, col) in izip!(trial_dofs, local_mat.col_iter()) { + for (test_dof, entry) in izip!(test_dofs, col.iter()) { + output.rows.push(test_space.global_dof_index(*test_dof)); + output.cols.push(trial_space.global_dof_index(*trial_dof)); + output.data.push(entry); + } + } + } + output +} + +fn get_pairs_if_smallest( + test_cell: &impl Entity, + trial_cell: &impl Entity, + vertex: usize, +) -> Option> { + let mut pairs = vec![]; + for (trial_i, trial_v) in trial_cell.topology().sub_entity_iter(0).enumerate() { + for (test_i, test_v) in test_cell.topology().sub_entity_iter(0).enumerate() { + if test_v == trial_v { + if test_v < vertex { + return None; + } + pairs.push((test_i, trial_i)); + } + } + } + Some(pairs) +} + +/// Options for a boundary assembler +pub struct BoundaryAssemblerOptions { + /// Number of points used in quadrature for non-singular integrals + quadrature_degrees: HashMap, + /// Quadrature degrees to be used for singular integrals + singular_quadrature_degrees: HashMap<(ReferenceCellType, ReferenceCellType), usize>, + /// Maximum size of each batch of cells to send to an assembly function + batch_size: usize, +} + +impl Default for BoundaryAssemblerOptions { + fn default() -> Self { + use ReferenceCellType::{Quadrilateral, Triangle}; + Self { + quadrature_degrees: HashMap::from([(Triangle, 37), (Quadrilateral, 37)]), + singular_quadrature_degrees: HashMap::from([ + ((Triangle, Triangle), 4), + ((Quadrilateral, Quadrilateral), 4), + ((Quadrilateral, Triangle), 4), + ((Triangle, Quadrilateral), 4), + ]), + batch_size: 128, + } + } +} + +// TODO: make this a struct, with HasBoundaryAssemblerOptions trait +pub trait BoundaryAssembler: Sync + Sized { + //! Boundary assembler + //! + //! Assemble operators by processing batches of cells in parallel + + /// Scalar type + type T: RlstScalar + MatrixInverse; + /// Integrand type + type Integrand: BoundaryIntegrand; + /// Kernel type + type Kernel: KernelEvaluator; + /// Number of derivatives + const DERIV_SIZE: usize; + /// Number of derivatives needed in basis function tables + const TABLE_DERIVS: usize; + + /// Get integrand + fn integrand(&self) -> &Self::Integrand; + + /// Get integrand + fn kernel(&self) -> &Self::Kernel; + + /// Get assembler options + fn options(&self) -> &BoundaryAssemblerOptions; + + /// Get mutable assembler options + fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions; + + /// Set (non-singular) quadrature degree for a cell type + fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) { + *self + .options_mut() + .quadrature_degrees + .get_mut(&cell) + .unwrap() = degree; + } + + /// Set singular quadrature degree for a pair of cell types + fn singular_quadrature_degree( + &mut self, + cells: (ReferenceCellType, ReferenceCellType), + degree: usize, + ) { + *self + .options_mut() + .singular_quadrature_degrees + .get_mut(&cells) + .unwrap() = degree; + } + + /// Set the maximum size of a batch of cells to send to an assembly function + fn batch_size(&mut self, size: usize) { + self.options_mut().batch_size = size; + } +} + +trait InternalAssemblyFunctions: BoundaryAssembler { + /// Assemble the singular contributions + fn assemble_singular< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + shape: [usize; 2], + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> SparseMatrixData { + if !equal_grids(test_space.grid(), trial_space.grid()) { + // If the test and trial grids are different, there are no neighbouring triangles + return SparseMatrixData::new(shape); + } + + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assembly can only be used for function spaces stored in serial"); + } + if shape[0] != test_space.global_size() || shape[1] != trial_space.global_size() { + panic!("Matrix has wrong shape"); + } + + let grid = test_space.grid(); + + let mut qweights = vec![]; + let mut trial_points = vec![]; + let mut test_points = vec![]; + let mut trial_tables = vec![]; + let mut test_tables = vec![]; + let mut test_cell_types = vec![]; + let mut trial_cell_types = vec![]; + + let mut pair_indices = HashMap::new(); + + for test_cell_type in grid.entity_types(2) { + for trial_cell_type in grid.entity_types(2) { + let qdegree = self.options().singular_quadrature_degrees + [&(*test_cell_type, *trial_cell_type)]; + let offset = qweights.len(); + + let mut possible_pairs = vec![]; + // Vertex-adjacent + for i in 0..reference_cell::entity_counts(*test_cell_type)[0] { + for j in 0..reference_cell::entity_counts(*trial_cell_type)[0] { + possible_pairs.push(vec![(i, j)]); + } + } + // edge-adjacent + for test_e in reference_cell::edges(*test_cell_type) { + for trial_e in reference_cell::edges(*trial_cell_type) { + possible_pairs.push(vec![(test_e[0], trial_e[0]), (test_e[1], trial_e[1])]); + possible_pairs.push(vec![(test_e[1], trial_e[0]), (test_e[0], trial_e[1])]); + } + } + // Same cell + if test_cell_type == trial_cell_type { + possible_pairs.push( + (0..reference_cell::entity_counts(*test_cell_type)[0]) + .map(&|i| (i, i)) + .collect::>(), + ); + } + + for (i, pairs) in possible_pairs.iter().enumerate() { + pair_indices.insert( + (*test_cell_type, *trial_cell_type, pairs.clone()), + offset + i, + ); + test_cell_types.push(*test_cell_type); + trial_cell_types.push(*trial_cell_type); + } + + for pairs in &possible_pairs { + let qrule = get_singular_quadrature_rule( + *test_cell_type, + *trial_cell_type, + pairs, + qdegree, + ); + let npts = qrule.weights.len(); + + let mut points = rlst_dynamic_array2!(::Real, [2, npts]); + for i in 0..npts { + for j in 0..2 { + *points.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule.trial_points[2 * i + j], + ) + .unwrap(); + } + } + let trial_element = trial_space.element(*trial_cell_type); + let mut table = rlst_dynamic_array4!( + Self::T, + trial_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1]) + ); + trial_element.tabulate(&points, Self::TABLE_DERIVS, &mut table); + trial_points.push(points); + trial_tables.push(table); + + let mut points = rlst_dynamic_array2!(::Real, [2, npts]); + for i in 0..npts { + for j in 0..2 { + *points.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule.test_points[2 * i + j], + ) + .unwrap(); + } + } + let test_element = test_space.element(*test_cell_type); + let mut table = rlst_dynamic_array4!( + Self::T, + test_element.tabulate_array_shape(Self::TABLE_DERIVS, points.shape()[1]) + ); + test_element.tabulate(&points, Self::TABLE_DERIVS, &mut table); + test_points.push(points); + test_tables.push(table); + qweights.push( + qrule + .weights + .iter() + .map(|w| num::cast::::Real>(*w).unwrap()) + .collect::>(), + ); + } + } + } + let cell_blocks = make_cell_blocks( + |test_cell_type, trial_cell_type, pairs| { + pair_indices[&(test_cell_type, trial_cell_type, pairs)] + }, + pair_indices.len(), + grid, + self.options().batch_size, + ); + + cell_blocks + .into_par_iter() + .map(|(i, cell_block)| { + assemble_batch_singular::( + self, + Self::DERIV_SIZE, + shape, + trial_cell_types[i], + test_cell_types[i], + trial_space, + test_space, + &cell_block, + &trial_points[i], + &test_points[i], + &qweights[i], + &trial_tables[i], + &test_tables[i], + ) + }) + .reduce( + || SparseMatrixData::::new(shape), + |mut a, b| { + a.add(b); + a + }, + ) + } + + /// Assemble the singular correction + /// + /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule + fn assemble_singular_correction< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + shape: [usize; 2], + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> SparseMatrixData { + if !equal_grids(test_space.grid(), trial_space.grid()) { + // If the test and trial grids are different, there are no neighbouring triangles + return SparseMatrixData::new(shape); + } + + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assembly can only be used for function spaces stored in serial"); + } + if shape[0] != test_space.global_size() || shape[1] != trial_space.global_size() { + panic!("Matrix has wrong shape"); + } + + let grid = test_space.grid(); + + let mut qweights_test = vec![]; + let mut qweights_trial = vec![]; + let mut qpoints_test = vec![]; + let mut qpoints_trial = vec![]; + let mut test_tables = vec![]; + let mut trial_tables = vec![]; + let mut test_cell_types = vec![]; + let mut trial_cell_types = vec![]; + + let mut cell_type_indices = HashMap::new(); + + for test_cell_type in grid.entity_types(2) { + let npts_test = self.options().quadrature_degrees[test_cell_type]; + for trial_cell_type in grid.entity_types(2) { + let npts_trial = self.options().quadrature_degrees[trial_cell_type]; + test_cell_types.push(*test_cell_type); + trial_cell_types.push(*trial_cell_type); + cell_type_indices.insert((*test_cell_type, *trial_cell_type), qweights_test.len()); + + let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap(); + let mut test_pts = + rlst_dynamic_array2!(::Real, [2, npts_test]); + for i in 0..npts_test { + for j in 0..2 { + *test_pts.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule_test.points[2 * i + j], + ) + .unwrap(); + } + } + qweights_test.push( + qrule_test + .weights + .iter() + .map(|w| num::cast::::Real>(*w).unwrap()) + .collect::>(), + ); + let test_element = test_space.element(*test_cell_type); + let mut test_table = rlst_dynamic_array4!( + Self::T, + test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test) + ); + test_element.tabulate(&test_pts, Self::TABLE_DERIVS, &mut test_table); + test_tables.push(test_table); + qpoints_test.push(test_pts); + + let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap(); + let mut trial_pts = + rlst_dynamic_array2!(::Real, [2, npts_trial]); + for i in 0..npts_trial { + for j in 0..2 { + *trial_pts.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule_trial.points[2 * i + j], + ) + .unwrap(); + } + } + qweights_trial.push( + qrule_trial + .weights + .iter() + .map(|w| num::cast::::Real>(*w).unwrap()) + .collect::>(), + ); + let trial_element = trial_space.element(*trial_cell_type); + let mut trial_table = rlst_dynamic_array4!( + Self::T, + trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial) + ); + trial_element.tabulate(&trial_pts, Self::TABLE_DERIVS, &mut trial_table); + trial_tables.push(trial_table); + qpoints_trial.push(trial_pts); + } + } + + let cell_blocks = make_cell_blocks( + |test_cell_type, trial_cell_type, pairs| { + cell_type_indices[&(test_cell_type, trial_cell_type)] + }, + qweights_test.len(), + grid, + self.options().batch_size, + ); + + cell_blocks + .into_par_iter() + .map(|(i, cell_block)| { + assemble_batch_singular_correction::( + self, + Self::DERIV_SIZE, + shape, + trial_cell_types[i], + test_cell_types[i], + trial_space, + test_space, + &cell_block, + &qpoints_trial[i], + &qweights_trial[i], + &qpoints_test[i], + &qweights_test[i], + &trial_tables[i], + &test_tables[i], + ) + }) + .reduce( + || SparseMatrixData::::new(shape), + |mut a, b| { + a.add(b); + a + }, + ) + } + /// Assemble the non-singular contributions into a dense matrix + fn assemble_nonsingular< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &RawData2D, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + trial_colouring: &HashMap>>, + test_colouring: &HashMap>>, + ) { + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assembly can only be used for function spaces stored in serial"); + } + if output.shape[0] != test_space.global_size() + || output.shape[1] != trial_space.global_size() + { + panic!("Matrix has wrong shape"); + } + + let batch_size = self.options().batch_size; + + for test_cell_type in test_space.grid().entity_types(2) { + let npts_test = self.options().quadrature_degrees[test_cell_type]; + for trial_cell_type in trial_space.grid().entity_types(2) { + let npts_trial = self.options().quadrature_degrees[trial_cell_type]; + let qrule_test = simplex_rule(*test_cell_type, npts_test).unwrap(); + let mut qpoints_test = + rlst_dynamic_array2!(::Real, [2, npts_test]); + for i in 0..npts_test { + for j in 0..2 { + *qpoints_test.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule_test.points[2 * i + j], + ) + .unwrap(); + } + } + let qweights_test = qrule_test + .weights + .iter() + .map(|w| num::cast::::Real>(*w).unwrap()) + .collect::>(); + let qrule_trial = simplex_rule(*trial_cell_type, npts_trial).unwrap(); + let mut qpoints_trial = + rlst_dynamic_array2!(::Real, [2, npts_trial]); + for i in 0..npts_trial { + for j in 0..2 { + *qpoints_trial.get_mut([j, i]).unwrap() = + num::cast::::Real>( + qrule_trial.points[2 * i + j], + ) + .unwrap(); + } + } + let qweights_trial = qrule_trial + .weights + .iter() + .map(|w| num::cast::::Real>(*w).unwrap()) + .collect::>(); + + let test_element = test_space.element(*test_cell_type); + let mut test_table = rlst_dynamic_array4!( + Self::T, + test_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_test) + ); + test_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut test_table); + + let trial_element = trial_space.element(*trial_cell_type); + let mut trial_table = rlst_dynamic_array4!( + Self::T, + trial_element.tabulate_array_shape(Self::TABLE_DERIVS, npts_trial) + ); + trial_element.tabulate(&qpoints_test, Self::TABLE_DERIVS, &mut trial_table); + + for test_c in &test_colouring[test_cell_type] { + for trial_c in &trial_colouring[trial_cell_type] { + let mut test_cells: Vec<&[usize]> = vec![]; + let mut trial_cells: Vec<&[usize]> = vec![]; + + let mut test_start = 0; + while test_start < test_c.len() { + let test_end = if test_start + batch_size < test_c.len() { + test_start + batch_size + } else { + test_c.len() + }; + + let mut trial_start = 0; + while trial_start < trial_c.len() { + let trial_end = if trial_start + batch_size < trial_c.len() { + trial_start + batch_size + } else { + trial_c.len() + }; + test_cells.push(&test_c[test_start..test_end]); + trial_cells.push(&trial_c[trial_start..trial_end]); + trial_start = trial_end; + } + test_start = test_end + } + + let numtasks = test_cells.len(); + let r: usize = (0..numtasks) + .into_par_iter() + .map(&|t| { + assemble_batch_nonadjacent::( + self, + Self::DERIV_SIZE, + output, + *test_cell_type, + *trial_cell_type, + trial_space, + trial_cells[t], + test_space, + test_cells[t], + &qpoints_trial, + &qweights_trial, + &qpoints_test, + &qweights_test, + &trial_table, + &test_table, + ) + }) + .sum(); + assert_eq!(r, numtasks); + } + } + } + } + } +} + +use crate::traits::BoundaryAssembly; +#[cfg(feature = "mpi")] +use crate::traits::ParallelBoundaryAssembly; + +impl InternalAssemblyFunctions for A {} +impl BoundaryAssembly for A { + type T = A::T; + fn assemble_singular_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) { + let sparse_matrix = self.assemble_singular::( + output.shape(), + trial_space, + test_space, + ); + let data = sparse_matrix.data; + let rows = sparse_matrix.rows; + let cols = sparse_matrix.cols; + for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { + *output.get_mut([*i, *j]).unwrap() += *value; + } + } + + fn assemble_singular_into_csr< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> CsrMatrix { + let shape = [test_space.global_size(), trial_space.global_size()]; + let sparse_matrix = + self.assemble_singular::(shape, trial_space, test_space); + + CsrMatrix::::from_aij( + sparse_matrix.shape, + &sparse_matrix.rows, + &sparse_matrix.cols, + &sparse_matrix.data, + ) + .unwrap() + } + + fn assemble_singular_correction_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) { + let sparse_matrix = self.assemble_singular_correction::( + output.shape(), + trial_space, + test_space, + ); + let data = sparse_matrix.data; + let rows = sparse_matrix.rows; + let cols = sparse_matrix.cols; + for ((i, j), value) in rows.iter().zip(cols.iter()).zip(data.iter()) { + *output.get_mut([*i, *j]).unwrap() += *value; + } + } + + fn assemble_singular_correction_into_csr< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> CsrMatrix { + let shape = [test_space.global_size(), trial_space.global_size()]; + let sparse_matrix = self.assemble_singular_correction::( + shape, + trial_space, + test_space, + ); + + CsrMatrix::::from_aij( + sparse_matrix.shape, + &sparse_matrix.rows, + &sparse_matrix.cols, + &sparse_matrix.data, + ) + .unwrap() + } + + fn assemble_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) { + let test_colouring = test_space.cell_colouring(); + let trial_colouring = trial_space.cell_colouring(); + + self.assemble_nonsingular_into_dense::( + output, + trial_space, + test_space, + &trial_colouring, + &test_colouring, + ); + self.assemble_singular_into_dense::( + output, + trial_space, + test_space, + ); + } + + fn assemble_nonsingular_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + trial_colouring: &HashMap>>, + test_colouring: &HashMap>>, + ) { + if !trial_space.is_serial() || !test_space.is_serial() { + panic!("Dense assembly can only be used for function spaces stored in serial"); + } + if output.shape()[0] != test_space.global_size() + || output.shape()[1] != trial_space.global_size() + { + panic!("Matrix has wrong shape"); + } + + let output_raw = RawData2D { + data: output.data_mut().as_mut_ptr(), + shape: output.shape(), + }; + + self.assemble_nonsingular( + &output_raw, + trial_space, + test_space, + trial_colouring, + test_colouring, + ) + } +} +#[cfg(feature = "mpi")] +impl ParallelBoundaryAssembly for A { + fn parallel_assemble_singular_into_csr< + 'a, + C: Communicator, + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + SerialTestSpace: FunctionSpace + Sync + 'a, + SerialTrialSpace: FunctionSpace + Sync + 'a, + >( + &self, + trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), + test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + ) -> CsrMatrix { + self.assemble_singular_into_csr(trial_space.local_space(), test_space.local_space()) + } + + fn parallel_assemble_singular_correction_into_csr< + 'a, + C: Communicator, + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + SerialTestSpace: FunctionSpace + Sync + 'a, + SerialTrialSpace: FunctionSpace + Sync + 'a, + >( + &self, + trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), + test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + ) -> CsrMatrix { + self.assemble_singular_correction_into_csr( + trial_space.local_space(), + test_space.local_space(), + ) + } +} diff --git a/src/assembly/boundary/adjoint_double_layer.rs b/src/assembly/boundary/assemblers/adjoint_double_layer.rs similarity index 100% rename from src/assembly/boundary/adjoint_double_layer.rs rename to src/assembly/boundary/assemblers/adjoint_double_layer.rs diff --git a/src/assembly/boundary/double_layer.rs b/src/assembly/boundary/assemblers/double_layer.rs similarity index 100% rename from src/assembly/boundary/double_layer.rs rename to src/assembly/boundary/assemblers/double_layer.rs diff --git a/src/assembly/boundary/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs similarity index 100% rename from src/assembly/boundary/hypersingular.rs rename to src/assembly/boundary/assemblers/hypersingular.rs diff --git a/src/assembly/boundary/single_layer.rs b/src/assembly/boundary/assemblers/single_layer.rs similarity index 100% rename from src/assembly/boundary/single_layer.rs rename to src/assembly/boundary/assemblers/single_layer.rs diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs index d4323069..793cee87 100644 --- a/src/assembly/boundary/cell_pair_assemblers.rs +++ b/src/assembly/boundary/cell_pair_assemblers.rs @@ -85,8 +85,9 @@ impl< I: BoundaryIntegrand, G: GeometryMap, K: KernelEvaluator, - > CellPairAssembler for SingularCellPairAssembler<'a, T, I, G, K> + > CellPairAssembler for SingularCellPairAssembler<'a, T, I, G, K> { + type T = T; fn set_test_cell(&mut self, test_cell: usize) { self.test_cell = test_cell; self.test_evaluator @@ -239,8 +240,9 @@ impl< TestG: GeometryMap, TrialG: GeometryMap, K: KernelEvaluator, - > CellPairAssembler for NonsingularCellPairAssembler<'a, T, I, TestG, TrialG, K> + > CellPairAssembler for NonsingularCellPairAssembler<'a, T, I, TestG, TrialG, K> { + type T = T; fn set_test_cell(&mut self, test_cell: usize) { self.test_cell = test_cell; self.test_evaluator diff --git a/src/traits.rs b/src/traits.rs index 5e9c1b3e..5962fb4b 100644 --- a/src/traits.rs +++ b/src/traits.rs @@ -3,7 +3,11 @@ mod assembly; mod function; -pub use assembly::{BoundaryIntegrand, CellGeometry, CellPairAssembler, KernelEvaluator}; +#[cfg(feature = "mpi")] +pub use assembly::ParallelBoundaryAssembly; +pub use assembly::{ + BoundaryAssembly, BoundaryIntegrand, CellGeometry, CellPairAssembler, KernelEvaluator, +}; pub use function::FunctionSpace; #[cfg(feature = "mpi")] pub use function::ParallelFunctionSpace; diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs index 2a4cdbb9..9d6d3859 100644 --- a/src/traits/assembly.rs +++ b/src/traits/assembly.rs @@ -1,6 +1,10 @@ //! Assembly use crate::assembly::common::RlstArray; -use rlst::RlstScalar; +use crate::traits::FunctionSpace; +use ndelement::{traits::FiniteElement, types::ReferenceCellType}; +use ndgrid::traits::Grid; +use rlst::{CsrMatrix, RlstScalar}; +use std::collections::HashMap; pub trait CellGeometry { //! Cell geometry @@ -75,12 +79,132 @@ pub trait KernelEvaluator { ); } -pub trait CellPairAssembler { +pub trait CellPairAssembler { //! Assembler for the contributions from a pair of cells + /// Scalar type + type T: RlstScalar; + /// Assemble contributions into `local_mat` - fn assemble(&mut self, local_mat: &mut RlstArray); + fn assemble(&mut self, local_mat: &mut RlstArray); /// Set the test cell fn set_test_cell(&mut self, test_cell: usize); /// Set the trial cell fn set_trial_cell(&mut self, trial_cell: usize); } + +pub trait BoundaryAssembly { + //! Functions for boundary assembly + /// Scalar type + type T: RlstScalar; + + /// Assemble the singular contributions into a dense matrix + fn assemble_singular_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ); + + /// Assemble the singular contributions into a CSR sparse matrix + fn assemble_singular_into_csr< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> CsrMatrix; + + /// Assemble the singular correction into a dense matrix + /// + /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule + fn assemble_singular_correction_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ); + + /// Assemble the singular correction into a CSR matrix + /// + /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule + fn assemble_singular_correction_into_csr< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ) -> CsrMatrix; + + /// Assemble into a dense matrix + fn assemble_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + ); + + /// Assemble the non-singular contributions into a dense matrix + fn assemble_nonsingular_into_dense< + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + >( + &self, + output: &mut RlstArray, + trial_space: &(impl FunctionSpace + Sync), + test_space: &(impl FunctionSpace + Sync), + trial_colouring: &HashMap>>, + test_colouring: &HashMap>>, + ); +} + +#[cfg(feature = "mpi")] +pub trait ParallelBoundaryAssembly: BoundaryAssembly { + //! Functions for parallel boundary assembly + + /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers + fn parallel_assemble_singular_into_csr< + 'a, + C: Communicator, + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + SerialTestSpace: FunctionSpace + Sync + 'a, + SerialTrialSpace: FunctionSpace + Sync + 'a, + >( + &self, + trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), + test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + ) -> CsrMatrix; + + /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers + fn parallel_assemble_singular_correction_into_csr< + 'a, + C: Communicator, + TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, + Element: FiniteElement + Sync, + SerialTestSpace: FunctionSpace + Sync + 'a, + SerialTrialSpace: FunctionSpace + Sync + 'a, + >( + &self, + trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), + test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + ) -> CsrMatrix; +} diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs index 308ec67d..adb4c132 100644 --- a/tests/compare_to_bempp_cl.rs +++ b/tests/compare_to_bempp_cl.rs @@ -3,7 +3,7 @@ use bempp::assembly::{ boundary, boundary::BoundaryAssembler, potential, potential::PotentialAssembler, }; use bempp::function::SerialFunctionSpace; -use bempp::traits::FunctionSpace; +use bempp::traits::{BoundaryAssembly, FunctionSpace}; use cauchy::c64; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::Continuity; From 1d911210609e0b938251182dd30df3f349bbd14a Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 13:53:25 +0100 Subject: [PATCH 06/13] simplify parameters --- src/assembly/boundary/assemblers.rs | 126 +++++++----------------- src/assembly/fmm_tools.rs | 12 +-- src/function/function_space/common.rs | 6 +- src/function/function_space/parallel.rs | 13 +-- src/function/function_space/serial.rs | 7 +- src/traits/assembly.rs | 60 ++++------- src/traits/function.rs | 7 +- 7 files changed, 80 insertions(+), 151 deletions(-) diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index d758bdb6..7624ebcf 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -499,15 +499,11 @@ pub trait BoundaryAssembler: Sync + Sized { trait InternalAssemblyFunctions: BoundaryAssembler { /// Assemble the singular contributions - fn assemble_singular< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular + Sync>( &self, shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> SparseMatrixData { if !equal_grids(test_space.grid(), trial_space.grid()) { // If the test and trial grids are different, there are no neighbouring triangles @@ -639,7 +635,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler { cell_blocks .into_par_iter() .map(|(i, cell_block)| { - assemble_batch_singular::( + assemble_batch_singular( self, Self::DERIV_SIZE, shape, @@ -667,15 +663,11 @@ trait InternalAssemblyFunctions: BoundaryAssembler { /// Assemble the singular correction /// /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_correction + Sync>( &self, shape: [usize; 2], - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> SparseMatrixData { if !equal_grids(test_space.grid(), trial_space.grid()) { // If the test and trial grids are different, there are no neighbouring triangles @@ -780,7 +772,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler { cell_blocks .into_par_iter() .map(|(i, cell_block)| { - assemble_batch_singular_correction::( + assemble_batch_singular_correction( self, Self::DERIV_SIZE, shape, @@ -806,15 +798,11 @@ trait InternalAssemblyFunctions: BoundaryAssembler { ) } /// Assemble the non-singular contributions into a dense matrix - fn assemble_nonsingular< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_nonsingular + Sync>( &self, output: &RawData2D, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, trial_colouring: &HashMap>>, test_colouring: &HashMap>>, ) { @@ -913,7 +901,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler { let r: usize = (0..numtasks) .into_par_iter() .map(&|t| { - assemble_batch_nonadjacent::( + assemble_batch_nonadjacent( self, Self::DERIV_SIZE, output, @@ -947,21 +935,13 @@ use crate::traits::ParallelBoundaryAssembly; impl InternalAssemblyFunctions for A {} impl BoundaryAssembly for A { type T = A::T; - fn assemble_singular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) { - let sparse_matrix = self.assemble_singular::( - output.shape(), - trial_space, - test_space, - ); + let sparse_matrix = self.assemble_singular(output.shape(), trial_space, test_space); let data = sparse_matrix.data; let rows = sparse_matrix.rows; let cols = sparse_matrix.cols; @@ -970,18 +950,13 @@ impl BoundaryAssembly for A { } } - fn assemble_singular_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_into_csr + Sync>( &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix { let shape = [test_space.global_size(), trial_space.global_size()]; - let sparse_matrix = - self.assemble_singular::(shape, trial_space, test_space); + let sparse_matrix = self.assemble_singular(shape, trial_space, test_space); CsrMatrix::::from_aij( sparse_matrix.shape, @@ -992,21 +967,14 @@ impl BoundaryAssembly for A { .unwrap() } - fn assemble_singular_correction_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_correction_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) { - let sparse_matrix = self.assemble_singular_correction::( - output.shape(), - trial_space, - test_space, - ); + let sparse_matrix = + self.assemble_singular_correction(output.shape(), trial_space, test_space); let data = sparse_matrix.data; let rows = sparse_matrix.rows; let cols = sparse_matrix.cols; @@ -1015,21 +983,13 @@ impl BoundaryAssembly for A { } } - fn assemble_singular_correction_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_correction_into_csr + Sync>( &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix { let shape = [test_space.global_size(), trial_space.global_size()]; - let sparse_matrix = self.assemble_singular_correction::( - shape, - trial_space, - test_space, - ); + let sparse_matrix = self.assemble_singular_correction(shape, trial_space, test_space); CsrMatrix::::from_aij( sparse_matrix.shape, @@ -1040,42 +1000,30 @@ impl BoundaryAssembly for A { .unwrap() } - fn assemble_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) { let test_colouring = test_space.cell_colouring(); let trial_colouring = trial_space.cell_colouring(); - self.assemble_nonsingular_into_dense::( + self.assemble_nonsingular_into_dense( output, trial_space, test_space, &trial_colouring, &test_colouring, ); - self.assemble_singular_into_dense::( - output, - trial_space, - test_space, - ); + self.assemble_singular_into_dense(output, trial_space, test_space); } - fn assemble_nonsingular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_nonsingular_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, trial_colouring: &HashMap>>, test_colouring: &HashMap>>, ) { diff --git a/src/assembly/fmm_tools.rs b/src/assembly/fmm_tools.rs index edda0574..387fb60c 100644 --- a/src/assembly/fmm_tools.rs +++ b/src/assembly/fmm_tools.rs @@ -15,7 +15,7 @@ use rlst::{ /// Generate an array of all the quadrature points pub fn get_all_quadrature_points< T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( npts: usize, grid: &G, @@ -56,7 +56,7 @@ pub fn get_all_quadrature_points< pub fn basis_to_quadrature_into_dense< const BLOCKSIZE: usize, T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( output: &mut Array, 2>, 2>, npts: usize, @@ -75,7 +75,7 @@ pub fn basis_to_quadrature_into_dense< pub fn basis_to_quadrature_into_csr< const BLOCKSIZE: usize, T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( npts: usize, space: &SerialFunctionSpace<'_, T, G>, @@ -102,7 +102,7 @@ pub fn basis_to_quadrature_into_csr< pub fn transpose_basis_to_quadrature_into_dense< const BLOCKSIZE: usize, T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( output: &mut Array, 2>, 2>, npts: usize, @@ -122,7 +122,7 @@ pub fn transpose_basis_to_quadrature_into_dense< pub fn transpose_basis_to_quadrature_into_csr< const BLOCKSIZE: usize, T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( npts: usize, space: &SerialFunctionSpace<'_, T, G>, @@ -148,7 +148,7 @@ pub fn transpose_basis_to_quadrature_into_csr< fn basis_to_quadrature< const BLOCKSIZE: usize, T: RlstScalar + MatrixInverse, - G: Grid, + G: Grid + Sync, >( shape: [usize; 2], npts: usize, diff --git a/src/function/function_space/common.rs b/src/function/function_space/common.rs index 7e4de897..d8b72455 100644 --- a/src/function/function_space/common.rs +++ b/src/function/function_space/common.rs @@ -15,7 +15,7 @@ type OwnerData = Vec<(usize, usize, usize, usize)>; pub(crate) fn assign_dofs< T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, >( rank: usize, grid: &GridImpl, @@ -103,7 +103,7 @@ mod test { use ndgrid::shapes::{screen_quadrilaterals, screen_triangles}; fn run_test( - grid: &impl Grid, + grid: &(impl Grid + Sync), degree: usize, continuity: Continuity, ) { @@ -131,7 +131,7 @@ mod test { } fn run_test_rt( - grid: &impl Grid, + grid: &(impl Grid + Sync), degree: usize, continuity: Continuity, ) { diff --git a/src/function/function_space/parallel.rs b/src/function/function_space/parallel.rs index 3363609e..8b09e831 100644 --- a/src/function/function_space/parallel.rs +++ b/src/function/function_space/parallel.rs @@ -21,7 +21,7 @@ use std::collections::HashMap; pub struct LocalFunctionSpace< 'a, T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, > { serial_space: SerialFunctionSpace<'a, T, GridImpl>, global_size: usize, @@ -32,9 +32,10 @@ pub struct LocalFunctionSpace< impl< 'a, T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, > FunctionSpace for LocalFunctionSpace<'a, T, GridImpl> { + type T = T; type Grid = GridImpl; type FiniteElement = CiarletElement; @@ -72,7 +73,7 @@ pub struct ParallelFunctionSpace< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid, + GridImpl: ParallelGrid + Grid + Sync, > { grid: &'a GridImpl, local_space: LocalFunctionSpace<'a, T, >::LocalGrid<'a>>, @@ -82,7 +83,7 @@ impl< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid, + GridImpl: ParallelGrid + Grid + Sync, > ParallelFunctionSpace<'a, C, T, GridImpl> { /// Create new function space @@ -228,7 +229,7 @@ impl< 'g, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid, + GridImpl: ParallelGrid + Grid + Sync, > ParallelFunctionSpaceTrait for ParallelFunctionSpace<'g, C, T, GridImpl> { type ParallelGrid = GridImpl; @@ -254,7 +255,7 @@ impl< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid, + GridImpl: ParallelGrid + Grid + Sync, > FunctionSpace for ParallelFunctionSpace<'a, C, T, GridImpl> { type Grid = GridImpl; diff --git a/src/function/function_space/serial.rs b/src/function/function_space/serial.rs index 500b7210..76041767 100644 --- a/src/function/function_space/serial.rs +++ b/src/function/function_space/serial.rs @@ -16,7 +16,7 @@ use std::collections::HashMap; pub struct SerialFunctionSpace< 'a, T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, > { pub(crate) grid: &'a GridImpl, pub(crate) elements: HashMap>, @@ -28,7 +28,7 @@ pub struct SerialFunctionSpace< impl< 'a, T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, > SerialFunctionSpace<'a, T, GridImpl> { /// Create new function space @@ -60,9 +60,10 @@ impl< impl< 'a, T: RlstScalar + MatrixInverse, - GridImpl: Grid, + GridImpl: Grid + Sync, > FunctionSpace for SerialFunctionSpace<'a, T, GridImpl> { + type T = T; type Grid = GridImpl; type FiniteElement = CiarletElement; diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs index 9d6d3859..044712b6 100644 --- a/src/traits/assembly.rs +++ b/src/traits/assembly.rs @@ -98,77 +98,53 @@ pub trait BoundaryAssembly { type T: RlstScalar; /// Assemble the singular contributions into a dense matrix - fn assemble_singular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ); /// Assemble the singular contributions into a CSR sparse matrix - fn assemble_singular_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_into_csr + Sync>( &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix; /// Assemble the singular correction into a dense matrix /// /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_correction_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ); /// Assemble the singular correction into a CSR matrix /// /// The singular correction is the contribution is the terms for adjacent cells are assembled using an (incorrect) non-singular quadrature rule - fn assemble_singular_correction_into_csr< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_singular_correction_into_csr + Sync>( &self, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix; /// Assemble into a dense matrix - fn assemble_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, ); /// Assemble the non-singular contributions into a dense matrix - fn assemble_nonsingular_into_dense< - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - >( + fn assemble_nonsingular_into_dense + Sync>( &self, output: &mut RlstArray, - trial_space: &(impl FunctionSpace + Sync), - test_space: &(impl FunctionSpace + Sync), + trial_space: &Space, + test_space: &Space, trial_colouring: &HashMap>>, test_colouring: &HashMap>>, ); diff --git a/src/traits/function.rs b/src/traits/function.rs index f5643e1a..bd21a278 100644 --- a/src/traits/function.rs +++ b/src/traits/function.rs @@ -5,14 +5,17 @@ use ndelement::{traits::FiniteElement, types::ReferenceCellType}; #[cfg(feature = "mpi")] use ndgrid::traits::ParallelGrid; use ndgrid::{traits::Grid, types::Ownership}; +use rlst::RlstScalar; use std::collections::HashMap; /// A function space pub trait FunctionSpace { + /// Scalar type + type T: RlstScalar; /// The grid type - type Grid: Grid; + type Grid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync; /// The finite element type - type FiniteElement: FiniteElement; + type FiniteElement: FiniteElement + Sync; /// Get the grid that the element is defined on fn grid(&self) -> &Self::Grid; From 396feeeadfc9443db9e6c8b3a50ce5d00abd9c67 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 15:13:37 +0100 Subject: [PATCH 07/13] working on hypersingular --- examples/test_parallel_assembly.rs | 1 - src/assembly.rs | 1 - src/assembly/boundary.rs | 2 - src/assembly/boundary/assemblers.rs | 31 ++-- .../boundary/assemblers/hypersingular.rs | 151 ++-------------- src/assembly/boundary/integrands.rs | 102 ++++++++++- .../boundary/integrands/hypersingular.rs | 169 +++++++++++++++++- src/assembly/common.rs | 13 -- src/function/function_space/parallel.rs | 11 +- src/traits/assembly.rs | 30 ++-- src/traits/function.rs | 7 +- tests/compare_to_bempp_cl.rs | 4 +- 12 files changed, 308 insertions(+), 214 deletions(-) diff --git a/examples/test_parallel_assembly.rs b/examples/test_parallel_assembly.rs index c2bea395..f92357b2 100644 --- a/examples/test_parallel_assembly.rs +++ b/examples/test_parallel_assembly.rs @@ -5,7 +5,6 @@ use approx::assert_relative_eq; #[cfg(feature = "mpi")] use bempp::{ assembly::boundary, - assembly::boundary::BoundaryAssembler, function::{ParallelFunctionSpace, SerialFunctionSpace}, traits::{BoundaryAssembly, FunctionSpace, ParallelBoundaryAssembly}, }; diff --git a/src/assembly.rs b/src/assembly.rs index 941d1cbf..5426a2ff 100644 --- a/src/assembly.rs +++ b/src/assembly.rs @@ -7,7 +7,6 @@ pub mod potential; #[cfg(test)] mod test { - use super::boundary::BoundaryAssembler; use super::*; use crate::function::SerialFunctionSpace; use crate::traits::{BoundaryAssembly, FunctionSpace}; diff --git a/src/assembly/boundary.rs b/src/assembly/boundary.rs index 922a56d5..ba7668d3 100644 --- a/src/assembly/boundary.rs +++ b/src/assembly/boundary.rs @@ -4,8 +4,6 @@ pub(crate) mod cell_pair_assemblers; mod integrands; pub use assemblers::BoundaryAssembler; -#[cfg(feature = "mpi")] -pub use assemblers::ParallelBoundaryAssembler; pub use assemblers::{ AdjointDoubleLayerAssembler, DoubleLayerAssembler, HypersingularAssembler, SingleLayerAssembler, }; diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 7624ebcf..9a1b9bc5 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -15,6 +15,9 @@ use crate::quadrature::duffy::{ }; use crate::quadrature::simplex_rules::simplex_rule; use crate::quadrature::types::{CellToCellConnectivity, TestTrialNumericalQuadratureDefinition}; +use crate::traits::BoundaryAssembly; +#[cfg(feature = "mpi")] +use crate::traits::ParallelBoundaryAssembly; #[cfg(feature = "mpi")] use crate::traits::ParallelFunctionSpace; use crate::traits::{BoundaryIntegrand, CellPairAssembler, FunctionSpace, KernelEvaluator}; @@ -761,7 +764,7 @@ trait InternalAssemblyFunctions: BoundaryAssembler { } let cell_blocks = make_cell_blocks( - |test_cell_type, trial_cell_type, pairs| { + |test_cell_type, trial_cell_type, _pairs| { cell_type_indices[&(test_cell_type, trial_cell_type)] }, qweights_test.len(), @@ -928,10 +931,6 @@ trait InternalAssemblyFunctions: BoundaryAssembler { } } -use crate::traits::BoundaryAssembly; -#[cfg(feature = "mpi")] -use crate::traits::ParallelBoundaryAssembly; - impl InternalAssemblyFunctions for A {} impl BoundaryAssembly for A { type T = A::T; @@ -1053,33 +1052,23 @@ impl BoundaryAssembly for A { #[cfg(feature = "mpi")] impl ParallelBoundaryAssembly for A { fn parallel_assemble_singular_into_csr< - 'a, C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, + Space: ParallelFunctionSpace, >( &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix { self.assemble_singular_into_csr(trial_space.local_space(), test_space.local_space()) } fn parallel_assemble_singular_correction_into_csr< - 'a, C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, + Space: ParallelFunctionSpace, >( &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix { self.assemble_singular_correction_into_csr( trial_space.local_space(), diff --git a/src/assembly/boundary/assemblers/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs index 0650e2d8..dd923cdf 100644 --- a/src/assembly/boundary/assemblers/hypersingular.rs +++ b/src/assembly/boundary/assemblers/hypersingular.rs @@ -1,70 +1,20 @@ //! Hypersingular assemblers use super::{BoundaryAssembler, BoundaryAssemblerOptions}; use crate::assembly::{ - boundary::integrands::HypersingularBoundaryIntegrand, common::GreenKernelEvalType, + boundary::integrands::{ + BoundaryIntegrandSum, HypersingularCurlCurlBoundaryIntegrand, + HypersingularNormalNormalBoundaryIntegrand, + }, + common::GreenKernelEvalType, kernels::KernelEvaluator, }; use green_kernels::{helmholtz_3d::Helmholtz3dKernel, laplace_3d::Laplace3dKernel, traits::Kernel}; use rlst::{MatrixInverse, RlstScalar}; -/* -#[allow(clippy::too_many_arguments)] -unsafe fn hyp_test_trial_product( - test_table: &RlstArray, - trial_table: &RlstArray, - test_jacobians: &RlstArray, - trial_jacobians: &RlstArray, - test_jdets: &[T::Real], - trial_jdets: &[T::Real], - test_point_index: usize, - trial_point_index: usize, - test_basis_index: usize, - trial_basis_index: usize, -) -> T { - let test0 = *test_table.get_unchecked([1, test_point_index, test_basis_index, 0]); - let test1 = *test_table.get_unchecked([2, test_point_index, test_basis_index, 0]); - let trial0 = *trial_table.get_unchecked([1, trial_point_index, trial_basis_index, 0]); - let trial1 = *trial_table.get_unchecked([2, trial_point_index, trial_basis_index, 0]); - - ((num::cast::(*test_jacobians.get_unchecked([3, test_point_index])).unwrap() - * test0 - - num::cast::(*test_jacobians.get_unchecked([0, test_point_index])).unwrap() - * test1) - * (num::cast::(*trial_jacobians.get_unchecked([3, trial_point_index])) - .unwrap() - * trial0 - - num::cast::(*trial_jacobians.get_unchecked([0, trial_point_index])) - .unwrap() - * trial1) - + (num::cast::(*test_jacobians.get_unchecked([4, test_point_index])).unwrap() - * test0 - - num::cast::(*test_jacobians.get_unchecked([1, test_point_index])) - .unwrap() - * test1) - * (num::cast::(*trial_jacobians.get_unchecked([4, trial_point_index])) - .unwrap() - * trial0 - - num::cast::(*trial_jacobians.get_unchecked([1, trial_point_index])) - .unwrap() - * trial1) - + (num::cast::(*test_jacobians.get_unchecked([5, test_point_index])).unwrap() - * test0 - - num::cast::(*test_jacobians.get_unchecked([2, test_point_index])) - .unwrap() - * test1) - * (num::cast::(*trial_jacobians.get_unchecked([5, trial_point_index])) - .unwrap() - * trial0 - - num::cast::(*trial_jacobians.get_unchecked([2, trial_point_index])) - .unwrap() - * trial1)) - / num::cast::(test_jdets[test_point_index] * trial_jdets[trial_point_index]) - .unwrap() -} -*/ /// Assembler for a hypersingular operator pub struct HypersingularAssembler> { kernel: KernelEvaluator, + integrand: HypersingularCurlCurlBoundaryIntegrand, options: BoundaryAssemblerOptions, } impl> HypersingularAssembler { @@ -72,6 +22,7 @@ impl> HypersingularAssembler) -> Self { Self { kernel, + integrand: HypersingularCurlCurlBoundaryIntegrand::new(), options: BoundaryAssemblerOptions::default(), } } @@ -97,13 +48,13 @@ impl + MatrixInverse> HypersingularAssembler BoundaryAssembler for HypersingularAssembler> { - const DERIV_SIZE: usize = 1; + const DERIV_SIZE: usize = 4; const TABLE_DERIVS: usize = 1; type T = T; - type Integrand = HypersingularBoundaryIntegrand; + type Integrand = HypersingularCurlCurlBoundaryIntegrand; type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularBoundaryIntegrand { - panic!(); + fn integrand(&self) -> &HypersingularCurlCurlBoundaryIntegrand { + &self.integrand } fn kernel(&self) -> &KernelEvaluator> { &self.kernel @@ -116,91 +67,15 @@ impl BoundaryAssembler } } -/// Assembler for curl-curl term of Helmholtz hypersingular operator -struct HelmholtzHypersingularCurlCurlAssembler<'a, T: RlstScalar + MatrixInverse> { - kernel: &'a KernelEvaluator>, - options: &'a BoundaryAssemblerOptions, -} -impl<'a, T: RlstScalar + MatrixInverse> - HelmholtzHypersingularCurlCurlAssembler<'a, T> -{ - /// Create a new assembler - pub fn new( - kernel: &'a KernelEvaluator>, - options: &'a BoundaryAssemblerOptions, - ) -> Self { - Self { kernel, options } - } -} -impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler - for HelmholtzHypersingularCurlCurlAssembler<'a, T> -{ - const DERIV_SIZE: usize = 1; - const TABLE_DERIVS: usize = 1; - type T = T; - type Integrand = HypersingularBoundaryIntegrand; - type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularBoundaryIntegrand { - panic!(); - } - fn kernel(&self) -> &KernelEvaluator> { - panic!(); - } - fn options(&self) -> &BoundaryAssemblerOptions { - self.options - } - fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { - panic!("Cannot get mutable options") - } -} - -/// Assembler for normal normal term of Helmholtz hypersingular boundary operator -struct HelmholtzHypersingularNormalNormalAssembler<'a, T: RlstScalar + MatrixInverse> { - kernel: &'a KernelEvaluator>, - options: &'a BoundaryAssemblerOptions, -} -impl<'a, T: RlstScalar + MatrixInverse> - HelmholtzHypersingularNormalNormalAssembler<'a, T> -{ - /// Create a new assembler - pub fn new( - kernel: &'a KernelEvaluator>, - options: &'a BoundaryAssemblerOptions, - ) -> Self { - Self { kernel, options } - } -} -impl<'a, T: RlstScalar + MatrixInverse> BoundaryAssembler - for HelmholtzHypersingularNormalNormalAssembler<'a, T> -{ - const DERIV_SIZE: usize = 1; - const TABLE_DERIVS: usize = 0; - type T = T; - type Integrand = HypersingularBoundaryIntegrand; - type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularBoundaryIntegrand { - panic!(); - } - fn kernel(&self) -> &KernelEvaluator> { - panic!(); - } - fn options(&self) -> &BoundaryAssemblerOptions { - self.options - } - fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { - panic!("Cannot get mutable options") - } -} - impl + MatrixInverse> BoundaryAssembler for HypersingularAssembler> { const DERIV_SIZE: usize = 1; const TABLE_DERIVS: usize = 1; type T = T; - type Integrand = HypersingularBoundaryIntegrand; + type Integrand = HypersingularCurlCurlBoundaryIntegrand; type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularBoundaryIntegrand { + fn integrand(&self) -> &HypersingularCurlCurlBoundaryIntegrand { panic!(); } fn kernel(&self) -> &KernelEvaluator> { diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs index c7751b6d..41c6afef 100644 --- a/src/assembly/boundary/integrands.rs +++ b/src/assembly/boundary/integrands.rs @@ -6,5 +6,105 @@ mod single_layer; pub use adjoint_double_layer::AdjointDoubleLayerBoundaryIntegrand; pub use double_layer::DoubleLayerBoundaryIntegrand; -pub use hypersingular::HypersingularBoundaryIntegrand; +pub use hypersingular::{ + HypersingularCurlCurlBoundaryIntegrand, HypersingularNormalNormalBoundaryIntegrand, +}; pub use single_layer::SingleLayerBoundaryIntegrand; + +use crate::assembly::common::RlstArray; +use crate::traits::{BoundaryIntegrand, CellGeometry}; +use rlst::RlstScalar; + +/// The sum of two integrands +pub struct BoundaryIntegrandSum< + T: RlstScalar, + I0: BoundaryIntegrand, + I1: BoundaryIntegrand, +> { + integrand0: I0, + integrand1: I1, +} + +impl, I1: BoundaryIntegrand> + BoundaryIntegrandSum +{ + pub fn new(integrand0: I0, integrand1: I1) -> Self { + Self { + integrand0, + integrand1, + } + } +} + +impl, I1: BoundaryIntegrand> BoundaryIntegrand + for BoundaryIntegrandSum +{ + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + self.integrand0.evaluate_nonsingular( + test_table, + trial_table, + test_point_index, + trial_point_index, + test_basis_index, + trial_basis_index, + k, + test_geometry, + trial_geometry, + ) * self.integrand1.evaluate_nonsingular( + test_table, + trial_table, + test_point_index, + trial_point_index, + test_basis_index, + trial_basis_index, + k, + test_geometry, + trial_geometry, + ) + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + self.integrand0.evaluate_singular( + test_table, + trial_table, + point_index, + test_basis_index, + trial_basis_index, + k, + test_geometry, + trial_geometry, + ) * self.integrand1.evaluate_singular( + test_table, + trial_table, + point_index, + test_basis_index, + trial_basis_index, + k, + test_geometry, + trial_geometry, + ) + } +} diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs index 933b976d..2699a1eb 100644 --- a/src/assembly/boundary/integrands/hypersingular.rs +++ b/src/assembly/boundary/integrands/hypersingular.rs @@ -1,13 +1,67 @@ //! Hypersingular integrand use crate::assembly::common::RlstArray; use crate::traits::{BoundaryIntegrand, CellGeometry}; -use rlst::RlstScalar; +use rlst::{RlstScalar, UnsafeRandomAccessByRef}; -pub struct HypersingularBoundaryIntegrand { +#[allow(clippy::too_many_arguments)] +unsafe fn hyp_test_trial_product( + test_table: &RlstArray, + trial_table: &RlstArray, + test_jacobians: &RlstArray, + trial_jacobians: &RlstArray, + test_jdets: &[T::Real], + trial_jdets: &[T::Real], + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, +) -> T { + let test0 = *test_table.get_unchecked([1, test_point_index, test_basis_index, 0]); + let test1 = *test_table.get_unchecked([2, test_point_index, test_basis_index, 0]); + let trial0 = *trial_table.get_unchecked([1, trial_point_index, trial_basis_index, 0]); + let trial1 = *trial_table.get_unchecked([2, trial_point_index, trial_basis_index, 0]); + + ((num::cast::(*test_jacobians.get_unchecked([3, test_point_index])).unwrap() + * test0 + - num::cast::(*test_jacobians.get_unchecked([0, test_point_index])).unwrap() + * test1) + * (num::cast::(*trial_jacobians.get_unchecked([3, trial_point_index])) + .unwrap() + * trial0 + - num::cast::(*trial_jacobians.get_unchecked([0, trial_point_index])) + .unwrap() + * trial1) + + (num::cast::(*test_jacobians.get_unchecked([4, test_point_index])).unwrap() + * test0 + - num::cast::(*test_jacobians.get_unchecked([1, test_point_index])) + .unwrap() + * test1) + * (num::cast::(*trial_jacobians.get_unchecked([4, trial_point_index])) + .unwrap() + * trial0 + - num::cast::(*trial_jacobians.get_unchecked([1, trial_point_index])) + .unwrap() + * trial1) + + (num::cast::(*test_jacobians.get_unchecked([5, test_point_index])).unwrap() + * test0 + - num::cast::(*test_jacobians.get_unchecked([2, test_point_index])) + .unwrap() + * test1) + * (num::cast::(*trial_jacobians.get_unchecked([5, trial_point_index])) + .unwrap() + * trial0 + - num::cast::(*trial_jacobians.get_unchecked([2, trial_point_index])) + .unwrap() + * trial1)) + / num::cast::(test_jdets[test_point_index] * trial_jdets[trial_point_index]) + .unwrap() +} + +pub struct HypersingularCurlCurlBoundaryIntegrand { _t: std::marker::PhantomData, } -impl HypersingularBoundaryIntegrand { +impl HypersingularCurlCurlBoundaryIntegrand { pub fn new() -> Self { Self { _t: std::marker::PhantomData, @@ -15,7 +69,76 @@ impl HypersingularBoundaryIntegrand { } } -impl BoundaryIntegrand for HypersingularBoundaryIntegrand { +impl BoundaryIntegrand for HypersingularCurlCurlBoundaryIntegrand { + type T = T; + + unsafe fn evaluate_nonsingular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + test_point_index: usize, + trial_point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + hyp_test_trial_product( + test_table, + trial_table, + test_geometry.jacobians(), + trial_geometry.jacobians(), + test_geometry.jdets(), + trial_geometry.jdets(), + test_point_index, + trial_point_index, + test_basis_index, + trial_basis_index, + ) * *k.get_unchecked([0, test_point_index, trial_point_index]) + } + + unsafe fn evaluate_singular( + &self, + test_table: &RlstArray, + trial_table: &RlstArray, + point_index: usize, + test_basis_index: usize, + trial_basis_index: usize, + k: &RlstArray, + test_geometry: &impl CellGeometry, + trial_geometry: &impl CellGeometry, + ) -> T { + hyp_test_trial_product( + test_table, + trial_table, + test_geometry.jacobians(), + trial_geometry.jacobians(), + test_geometry.jdets(), + trial_geometry.jdets(), + point_index, + point_index, + test_basis_index, + trial_basis_index, + ) * *k.get_unchecked([0, point_index]) + } +} + +pub struct HypersingularNormalNormalBoundaryIntegrand { + wavenumber: T::Real, + _t: std::marker::PhantomData, +} + +impl HypersingularNormalNormalBoundaryIntegrand { + pub fn new(wavenumber: T::Real) -> Self { + Self { + wavenumber, + _t: std::marker::PhantomData, + } + } +} + +impl BoundaryIntegrand for HypersingularNormalNormalBoundaryIntegrand { type T = T; unsafe fn evaluate_nonsingular( @@ -27,10 +150,28 @@ impl BoundaryIntegrand for HypersingularBoundaryIntegrand { test_basis_index: usize, trial_basis_index: usize, k: &RlstArray, - _test_geometry: &impl CellGeometry, + test_geometry: &impl CellGeometry, trial_geometry: &impl CellGeometry, ) -> T { - panic!(); + -*k.get_unchecked([0, test_point_index, trial_point_index]) + * num::cast::( + self.wavenumber.powi(2) + * (*trial_geometry + .normals() + .get_unchecked([0, trial_point_index]) + * *test_geometry.normals().get_unchecked([0, test_point_index]) + + *trial_geometry + .normals() + .get_unchecked([1, trial_point_index]) + * *test_geometry.normals().get_unchecked([1, test_point_index]) + + *trial_geometry + .normals() + .get_unchecked([2, trial_point_index]) + * *test_geometry.normals().get_unchecked([2, test_point_index])), + ) + .unwrap() + * *test_table.get_unchecked([0, test_point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, trial_point_index, trial_basis_index, 0]) } unsafe fn evaluate_singular( @@ -41,9 +182,21 @@ impl BoundaryIntegrand for HypersingularBoundaryIntegrand { test_basis_index: usize, trial_basis_index: usize, k: &RlstArray, - _test_geometry: &impl CellGeometry, + test_geometry: &impl CellGeometry, trial_geometry: &impl CellGeometry, ) -> T { - panic!(); + -*k.get_unchecked([0, point_index]) + * num::cast::( + self.wavenumber.powi(2) + * (*trial_geometry.normals().get_unchecked([0, point_index]) + * *test_geometry.normals().get_unchecked([0, point_index]) + + *trial_geometry.normals().get_unchecked([1, point_index]) + * *test_geometry.normals().get_unchecked([1, point_index]) + + *trial_geometry.normals().get_unchecked([2, point_index]) + * *test_geometry.normals().get_unchecked([2, point_index])), + ) + .unwrap() + * *test_table.get_unchecked([0, point_index, test_basis_index, 0]) + * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) } } diff --git a/src/assembly/common.rs b/src/assembly/common.rs index fe184e06..7e33da5a 100644 --- a/src/assembly/common.rs +++ b/src/assembly/common.rs @@ -63,19 +63,6 @@ impl SparseMatrixData { self.cols.extend(&other.cols); self.data.extend(&other.data); } - /// Compute the sum of this sparse matrix and another sparse matrix - pub fn sum(&self, other: SparseMatrixData) -> SparseMatrixData { - debug_assert!(self.shape[0] == other.shape[0]); - debug_assert!(self.shape[1] == other.shape[1]); - let mut out = SparseMatrixData::::new(self.shape); - out.rows.extend(&self.rows); - out.cols.extend(&self.cols); - out.data.extend(&self.data); - out.rows.extend(&other.rows); - out.cols.extend(&other.cols); - out.data.extend(&other.data); - out - } } unsafe impl Sync for SparseMatrixData {} diff --git a/src/function/function_space/parallel.rs b/src/function/function_space/parallel.rs index 8b09e831..f6bd11ce 100644 --- a/src/function/function_space/parallel.rs +++ b/src/function/function_space/parallel.rs @@ -73,7 +73,7 @@ pub struct ParallelFunctionSpace< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid + Sync, + GridImpl: ParallelGrid + Grid, > { grid: &'a GridImpl, local_space: LocalFunctionSpace<'a, T, >::LocalGrid<'a>>, @@ -83,7 +83,7 @@ impl< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid + Sync, + GridImpl: ParallelGrid + Grid, > ParallelFunctionSpace<'a, C, T, GridImpl> { /// Create new function space @@ -101,7 +101,7 @@ impl< // Create local space on current process let (cell_dofs, entity_dofs, dofmap_size, owner_data) = - assign_dofs(rank as usize, grid, e_family); + assign_dofs(rank as usize, grid.local_grid(), e_family); let mut elements = HashMap::new(); for cell in grid.entity_types(grid.topology_dim()) { @@ -229,7 +229,7 @@ impl< 'g, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid + Sync, + GridImpl: ParallelGrid + Grid, > ParallelFunctionSpaceTrait for ParallelFunctionSpace<'g, C, T, GridImpl> { type ParallelGrid = GridImpl; @@ -255,9 +255,10 @@ impl< 'a, C: Communicator, T: RlstScalar + MatrixInverse, - GridImpl: ParallelGrid + Grid + Sync, + GridImpl: ParallelGrid + Grid, > FunctionSpace for ParallelFunctionSpace<'a, C, T, GridImpl> { + type T = T; type Grid = GridImpl; type FiniteElement = CiarletElement; diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs index 044712b6..2dfdfee7 100644 --- a/src/traits/assembly.rs +++ b/src/traits/assembly.rs @@ -1,11 +1,15 @@ //! Assembly use crate::assembly::common::RlstArray; use crate::traits::FunctionSpace; -use ndelement::{traits::FiniteElement, types::ReferenceCellType}; -use ndgrid::traits::Grid; +use ndelement::types::ReferenceCellType; use rlst::{CsrMatrix, RlstScalar}; use std::collections::HashMap; +#[cfg(feature = "mpi")] +use crate::traits::ParallelFunctionSpace; +#[cfg(feature = "mpi")] +use mpi::traits::Communicator; + pub trait CellGeometry { //! Cell geometry /// Scalar type @@ -156,31 +160,21 @@ pub trait ParallelBoundaryAssembly: BoundaryAssembly { /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers fn parallel_assemble_singular_into_csr< - 'a, C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, + Space: ParallelFunctionSpace, >( &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix; /// Assemble the singular contributions into a CSR sparse matrix, indexed by global DOF numbers fn parallel_assemble_singular_correction_into_csr< - 'a, C: Communicator, - TestGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - TrialGrid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync, - Element: FiniteElement + Sync, - SerialTestSpace: FunctionSpace + Sync + 'a, - SerialTrialSpace: FunctionSpace + Sync + 'a, + Space: ParallelFunctionSpace, >( &self, - trial_space: &'a (impl ParallelFunctionSpace = SerialTrialSpace> + 'a), - test_space: &'a (impl ParallelFunctionSpace = SerialTestSpace> + 'a), + trial_space: &Space, + test_space: &Space, ) -> CsrMatrix; } diff --git a/src/traits/function.rs b/src/traits/function.rs index bd21a278..56a9f3ae 100644 --- a/src/traits/function.rs +++ b/src/traits/function.rs @@ -13,7 +13,7 @@ pub trait FunctionSpace { /// Scalar type type T: RlstScalar; /// The grid type - type Grid: Grid::Real, EntityDescriptor = ReferenceCellType> + Sync; + type Grid: Grid::Real, EntityDescriptor = ReferenceCellType>; /// The finite element type type FiniteElement: FiniteElement + Sync; @@ -55,9 +55,10 @@ pub trait FunctionSpace { /// A function space in parallel pub trait ParallelFunctionSpace: FunctionSpace { /// Parallel grid type - type ParallelGrid: ParallelGrid + Grid; + type ParallelGrid: ParallelGrid + + Grid::Real, EntityDescriptor = ReferenceCellType>; /// The type of the serial space on each process - type LocalSpace<'a>: FunctionSpace + type LocalSpace<'a>: FunctionSpace + Sync where Self: 'a; diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs index adb4c132..fcaa789f 100644 --- a/tests/compare_to_bempp_cl.rs +++ b/tests/compare_to_bempp_cl.rs @@ -1,7 +1,5 @@ use approx::*; -use bempp::assembly::{ - boundary, boundary::BoundaryAssembler, potential, potential::PotentialAssembler, -}; +use bempp::assembly::{boundary, potential, potential::PotentialAssembler}; use bempp::function::SerialFunctionSpace; use bempp::traits::{BoundaryAssembly, FunctionSpace}; use cauchy::c64; From d69fd1742df2988f07c4b84395e73f9380eabcba Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 15:26:41 +0100 Subject: [PATCH 08/13] hypersingular --- src/assembly.rs | 10 +++ .../boundary/assemblers/hypersingular.rs | 84 +++++++++---------- tests/compare_to_bempp_cl.rs | 6 +- 3 files changed, 55 insertions(+), 45 deletions(-) diff --git a/src/assembly.rs b/src/assembly.rs index 5426a2ff..7570c7ad 100644 --- a/src/assembly.rs +++ b/src/assembly.rs @@ -121,6 +121,16 @@ mod test { //}; } 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::[<$operator Assembler>]::<[<$dtype>], _>::new_laplace() diff --git a/src/assembly/boundary/assemblers/hypersingular.rs b/src/assembly/boundary/assemblers/hypersingular.rs index dd923cdf..e0fb5a1a 100644 --- a/src/assembly/boundary/assemblers/hypersingular.rs +++ b/src/assembly/boundary/assemblers/hypersingular.rs @@ -8,55 +8,77 @@ use crate::assembly::{ 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, + HypersingularNormalNormalBoundaryIntegrand, +>; + /// Assembler for a hypersingular operator -pub struct HypersingularAssembler> { +pub struct HypersingularAssembler< + T: RlstScalar + MatrixInverse, + K: Kernel, + I: BoundaryIntegrand + Sync, +> { kernel: KernelEvaluator, - integrand: HypersingularCurlCurlBoundaryIntegrand, + integrand: I, options: BoundaryAssemblerOptions, } -impl> HypersingularAssembler { + +impl, I: BoundaryIntegrand + Sync> + HypersingularAssembler +{ /// Create a new hypersingular assembler - pub fn new(kernel: KernelEvaluator) -> Self { + pub fn new(kernel: KernelEvaluator, integrand: I) -> Self { Self { kernel, - integrand: HypersingularCurlCurlBoundaryIntegrand::new(), + integrand, options: BoundaryAssemblerOptions::default(), } } } -impl HypersingularAssembler> { +impl + HypersingularAssembler, HypersingularCurlCurlBoundaryIntegrand> +{ /// Create a new Laplace hypersingular assembler pub fn new_laplace() -> Self { - Self::new(KernelEvaluator::new_laplace( - GreenKernelEvalType::ValueDeriv, - )) + Self::new( + KernelEvaluator::new_laplace(GreenKernelEvalType::ValueDeriv), + HypersingularCurlCurlBoundaryIntegrand::new(), + ) } } -impl + MatrixInverse> HypersingularAssembler> { +impl + MatrixInverse> + HypersingularAssembler, HelmholtzIntegrand> +{ /// Create a new Helmholtz hypersingular assembler pub fn new_helmholtz(wavenumber: T::Real) -> Self { - Self::new(KernelEvaluator::new_helmholtz( - wavenumber, - GreenKernelEvalType::ValueDeriv, - )) + Self::new( + KernelEvaluator::new_helmholtz(wavenumber, GreenKernelEvalType::ValueDeriv), + BoundaryIntegrandSum::new( + HypersingularCurlCurlBoundaryIntegrand::new(), + HypersingularNormalNormalBoundaryIntegrand::new(wavenumber), + ), + ) } } -impl BoundaryAssembler - for HypersingularAssembler> +impl, I: BoundaryIntegrand + Sync> + BoundaryAssembler for HypersingularAssembler { const DERIV_SIZE: usize = 4; const TABLE_DERIVS: usize = 1; type T = T; - type Integrand = HypersingularCurlCurlBoundaryIntegrand; - type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularCurlCurlBoundaryIntegrand { + type Integrand = I; + type Kernel = KernelEvaluator; + fn integrand(&self) -> &I { &self.integrand } - fn kernel(&self) -> &KernelEvaluator> { + fn kernel(&self) -> &KernelEvaluator { &self.kernel } fn options(&self) -> &BoundaryAssemblerOptions { @@ -66,25 +88,3 @@ impl BoundaryAssembler &mut self.options } } - -impl + MatrixInverse> BoundaryAssembler - for HypersingularAssembler> -{ - const DERIV_SIZE: usize = 1; - const TABLE_DERIVS: usize = 1; - type T = T; - type Integrand = HypersingularCurlCurlBoundaryIntegrand; - type Kernel = KernelEvaluator>; - fn integrand(&self) -> &HypersingularCurlCurlBoundaryIntegrand { - panic!(); - } - fn kernel(&self) -> &KernelEvaluator> { - panic!(); - } - fn options(&self) -> &BoundaryAssemblerOptions { - &self.options - } - fn options_mut(&mut self) -> &mut BoundaryAssemblerOptions { - &mut self.options - } -} diff --git a/tests/compare_to_bempp_cl.rs b/tests/compare_to_bempp_cl.rs index fcaa789f..83f6b791 100644 --- a/tests/compare_to_bempp_cl.rs +++ b/tests/compare_to_bempp_cl.rs @@ -86,7 +86,7 @@ fn test_laplace_hypersingular_dp0_dp0() { let ndofs = space.global_size(); let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = boundary::HypersingularAssembler::::new_laplace(); + let a = boundary::HypersingularAssembler::::new_laplace(); a.assemble_into_dense(&mut matrix, &space, &space); for i in 0..ndofs { @@ -105,7 +105,7 @@ fn test_laplace_hypersingular_p1_p1() { let ndofs = space.global_size(); let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]); - let a = boundary::HypersingularAssembler::::new_laplace(); + let a = boundary::HypersingularAssembler::::new_laplace(); a.assemble_into_dense(&mut matrix, &space, &space); // Compare to result from bempp-cl @@ -202,7 +202,7 @@ fn test_helmholtz_hypersingular_p1_p1() { let ndofs = space.global_size(); let mut matrix = rlst_dynamic_array2!(c64, [ndofs, ndofs]); - let a = boundary::HypersingularAssembler::::new_helmholtz(3.0); + let a = boundary::HypersingularAssembler::::new_helmholtz(3.0); a.assemble_into_dense(&mut matrix, &space, &space); // Compare to result from bempp-cl From ce3336fbc2382ace93848175391a4d4b9962b454 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 15:46:16 +0100 Subject: [PATCH 09/13] fix some typos --- src/assembly/boundary/assemblers.rs | 4 ++-- src/assembly/boundary/integrands.rs | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 9a1b9bc5..94c474ba 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -291,7 +291,7 @@ fn assemble_batch_nonadjacent< ); for trial_cell in trial_cells { - a.set_test_cell(*trial_cell); + a.set_trial_cell(*trial_cell); let trial_dofs = trial_space.cell_dofs(*trial_cell).unwrap(); for test_cell in test_cells { if neighbours(test_grid, trial_grid, *test_cell, *trial_cell) { @@ -380,7 +380,7 @@ fn assemble_batch_singular_correction< for (test_cell, trial_cell) in cell_pairs { a.set_test_cell(*test_cell); - a.set_test_cell(*trial_cell); + a.set_trial_cell(*trial_cell); a.assemble(&mut local_mat); diff --git a/src/assembly/boundary/integrands.rs b/src/assembly/boundary/integrands.rs index 41c6afef..05494cde 100644 --- a/src/assembly/boundary/integrands.rs +++ b/src/assembly/boundary/integrands.rs @@ -63,7 +63,7 @@ impl, I1: BoundaryIntegrand> k, test_geometry, trial_geometry, - ) * self.integrand1.evaluate_nonsingular( + ) + self.integrand1.evaluate_nonsingular( test_table, trial_table, test_point_index, @@ -96,7 +96,7 @@ impl, I1: BoundaryIntegrand> k, test_geometry, trial_geometry, - ) * self.integrand1.evaluate_singular( + ) + self.integrand1.evaluate_singular( test_table, trial_table, point_index, From 7266352a7100116ff503eef0472d9e1a6a97deea Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 16:03:19 +0100 Subject: [PATCH 10/13] fix bench --- benches/assembly_benchmark.rs | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs index 1306d765..b8f74557 100644 --- a/benches/assembly_benchmark.rs +++ b/benches/assembly_benchmark.rs @@ -1,6 +1,6 @@ use bempp::assembly::{boundary, boundary::BoundaryAssembler}; use bempp::function::SerialFunctionSpace; -use bempp::traits::FunctionSpace; +use bempp::traits::{FunctionSpace, BoundaryAssembly}; use criterion::{criterion_group, criterion_main, Criterion}; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::{Continuity, ReferenceCellType}; From 00c49f75cb063876ac25f91a505b051544a4ea79 Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 16:04:28 +0100 Subject: [PATCH 11/13] remove old file --- src/assembly/batched.rs | 117 ---------------------------------------- 1 file changed, 117 deletions(-) delete mode 100644 src/assembly/batched.rs diff --git a/src/assembly/batched.rs b/src/assembly/batched.rs deleted file mode 100644 index 7cf2b4c7..00000000 --- a/src/assembly/batched.rs +++ /dev/null @@ -1,117 +0,0 @@ -//! Batched dense assembly -use crate::assembly::common::SparseMatrixData; -use green_kernels::types::EvalType as GreenKernelEvalType; -use rlst::{Array, BaseArray, VectorContainer}; - -mod boundary; -pub use boundary::adjoint_double_layer::{ - HelmholtzAdjointDoubleLayerAssembler, LaplaceAdjointDoubleLayerAssembler, -}; -pub use boundary::double_layer::{HelmholtzDoubleLayerAssembler, LaplaceDoubleLayerAssembler}; -pub use boundary::hypersingular::{HelmholtzHypersingularAssembler, LaplaceHypersingularAssembler}; -pub use boundary::single_layer::{HelmholtzSingleLayerAssembler, LaplaceSingleLayerAssembler}; -pub use boundary::{BatchedAssembler, BatchedAssemblerOptions}; - -mod potential; -pub use potential::double_layer::{ - HelmholtzDoubleLayerPotentialAssembler, LaplaceDoubleLayerPotentialAssembler, -}; -pub use potential::single_layer::{ - HelmholtzSingleLayerPotentialAssembler, LaplaceSingleLayerPotentialAssembler, -}; -pub use potential::{BatchedPotentialAssembler, BatchedPotentialAssemblerOptions}; - -type RlstArray = Array, DIM>, DIM>; - -#[cfg(test)] -mod test { - use super::*; - use crate::function::SerialFunctionSpace; - use crate::traits::FunctionSpace; - use approx::*; - use ndelement::ciarlet::LagrangeElementFamily; - use ndelement::types::Continuity; - use ndgrid::shapes::regular_sphere; - use rlst::rlst_dynamic_array2; - use rlst::RandomAccessByRef; - - #[test] - fn test_singular_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 assembler = LaplaceSingleLayerAssembler::::default(); - assembler.assemble_singular_into_dense(&mut matrix, &space, &space); - let csr = assembler.assemble_singular_into_csr(&space, &space); - - let indptr = csr.indptr(); - let indices = csr.indices(); - let data = csr.data(); - - let mut row = 0; - for (i, j) in indices.iter().enumerate() { - while i >= indptr[row + 1] { - row += 1; - } - assert_relative_eq!(*matrix.get([row, *j]).unwrap(), data[i], epsilon = 1e-8); - } - } - - #[test] - fn test_singular_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 assembler = LaplaceSingleLayerAssembler::::default(); - assembler.assemble_singular_into_dense(&mut matrix, &space, &space); - let csr = assembler.assemble_singular_into_csr(&space, &space); - - let indptr = csr.indptr(); - let indices = csr.indices(); - let data = csr.data(); - - let mut row = 0; - for (i, j) in indices.iter().enumerate() { - while i >= indptr[row + 1] { - row += 1; - } - assert_relative_eq!(*matrix.get([row, *j]).unwrap(), data[i], epsilon = 1e-8); - } - } - - #[test] - fn test_singular_dp0_p1() { - let grid = regular_sphere::(0); - let element0 = LagrangeElementFamily::::new(0, Continuity::Discontinuous); - let element1 = LagrangeElementFamily::::new(1, Continuity::Standard); - let space0 = SerialFunctionSpace::new(&grid, &element0); - let space1 = SerialFunctionSpace::new(&grid, &element1); - - let ndofs0 = space0.global_size(); - let ndofs1 = space1.global_size(); - - let mut matrix = rlst_dynamic_array2!(f64, [ndofs1, ndofs0]); - let assembler = LaplaceSingleLayerAssembler::::default(); - assembler.assemble_singular_into_dense(&mut matrix, &space0, &space1); - let csr = assembler.assemble_singular_into_csr(&space0, &space1); - let indptr = csr.indptr(); - let indices = csr.indices(); - let data = csr.data(); - - let mut row = 0; - for (i, j) in indices.iter().enumerate() { - while i >= indptr[row + 1] { - row += 1; - } - assert_relative_eq!(*matrix.get([row, *j]).unwrap(), data[i], epsilon = 1e-8); - } - } -} From bf4a63ca837aca38fbf8ed651d88e8d771f1b10a Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 16:08:26 +0100 Subject: [PATCH 12/13] clippy --- benches/assembly_benchmark.rs | 2 +- src/assembly/boundary/assemblers.rs | 8 ++++---- src/assembly/boundary/cell_pair_assemblers.rs | 2 ++ src/assembly/boundary/integrands/adjoint_double_layer.rs | 6 ++++++ src/assembly/boundary/integrands/double_layer.rs | 6 ++++++ src/assembly/boundary/integrands/hypersingular.rs | 6 ++++++ src/assembly/boundary/integrands/single_layer.rs | 6 ++++++ src/traits/assembly.rs | 8 ++++++++ 8 files changed, 39 insertions(+), 5 deletions(-) diff --git a/benches/assembly_benchmark.rs b/benches/assembly_benchmark.rs index b8f74557..c7adac70 100644 --- a/benches/assembly_benchmark.rs +++ b/benches/assembly_benchmark.rs @@ -1,6 +1,6 @@ use bempp::assembly::{boundary, boundary::BoundaryAssembler}; use bempp::function::SerialFunctionSpace; -use bempp::traits::{FunctionSpace, BoundaryAssembly}; +use bempp::traits::{BoundaryAssembly, FunctionSpace}; use criterion::{criterion_group, criterion_main, Criterion}; use ndelement::ciarlet::LagrangeElementFamily; use ndelement::types::{Continuity, ReferenceCellType}; diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 94c474ba..4f2dae9b 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -278,8 +278,8 @@ fn assemble_batch_nonadjacent< trial_evaluator, test_table, trial_table, - &test_weights, - &trial_weights, + test_weights, + trial_weights, ); let mut local_mat = rlst_dynamic_array2!( @@ -366,8 +366,8 @@ fn assemble_batch_singular_correction< trial_evaluator, test_table, trial_table, - &test_weights, - &trial_weights, + test_weights, + trial_weights, ); let mut local_mat = rlst_dynamic_array2!( diff --git a/src/assembly/boundary/cell_pair_assemblers.rs b/src/assembly/boundary/cell_pair_assemblers.rs index 793cee87..45849b63 100644 --- a/src/assembly/boundary/cell_pair_assemblers.rs +++ b/src/assembly/boundary/cell_pair_assemblers.rs @@ -45,6 +45,7 @@ impl< K: KernelEvaluator, > SingularCellPairAssembler<'a, T, I, G, K> { + #[allow(clippy::too_many_arguments)] /// Create new pub fn new( npts: usize, @@ -195,6 +196,7 @@ impl< K: KernelEvaluator, > NonsingularCellPairAssembler<'a, T, I, TestG, TrialG, K> { + #[allow(clippy::too_many_arguments)] /// Create new pub fn new( npts_test: usize, diff --git a/src/assembly/boundary/integrands/adjoint_double_layer.rs b/src/assembly/boundary/integrands/adjoint_double_layer.rs index 049dbfe1..6561ec35 100644 --- a/src/assembly/boundary/integrands/adjoint_double_layer.rs +++ b/src/assembly/boundary/integrands/adjoint_double_layer.rs @@ -73,3 +73,9 @@ impl BoundaryIntegrand for AdjointDoubleLayerBoundaryIntegrand * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) } } + +impl Default for AdjointDoubleLayerBoundaryIntegrand { + fn default() -> Self { + Self::new() + } +} diff --git a/src/assembly/boundary/integrands/double_layer.rs b/src/assembly/boundary/integrands/double_layer.rs index 835b497e..3bfb4f82 100644 --- a/src/assembly/boundary/integrands/double_layer.rs +++ b/src/assembly/boundary/integrands/double_layer.rs @@ -83,3 +83,9 @@ impl BoundaryIntegrand for DoubleLayerBoundaryIntegrand { * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) } } + +impl Default for DoubleLayerBoundaryIntegrand { + fn default() -> Self { + Self::new() + } +} diff --git a/src/assembly/boundary/integrands/hypersingular.rs b/src/assembly/boundary/integrands/hypersingular.rs index 2699a1eb..9638700f 100644 --- a/src/assembly/boundary/integrands/hypersingular.rs +++ b/src/assembly/boundary/integrands/hypersingular.rs @@ -69,6 +69,12 @@ impl HypersingularCurlCurlBoundaryIntegrand { } } +impl Default for HypersingularCurlCurlBoundaryIntegrand { + fn default() -> Self { + Self::new() + } +} + impl BoundaryIntegrand for HypersingularCurlCurlBoundaryIntegrand { type T = T; diff --git a/src/assembly/boundary/integrands/single_layer.rs b/src/assembly/boundary/integrands/single_layer.rs index 2e9803c8..e78ff098 100644 --- a/src/assembly/boundary/integrands/single_layer.rs +++ b/src/assembly/boundary/integrands/single_layer.rs @@ -51,3 +51,9 @@ impl BoundaryIntegrand for SingleLayerBoundaryIntegrand { * *trial_table.get_unchecked([0, point_index, trial_basis_index, 0]) } } + +impl Default for SingleLayerBoundaryIntegrand { + fn default() -> Self { + Self::new() + } +} diff --git a/src/traits/assembly.rs b/src/traits/assembly.rs index 2dfdfee7..ca8ff204 100644 --- a/src/traits/assembly.rs +++ b/src/traits/assembly.rs @@ -29,7 +29,11 @@ pub trait BoundaryIntegrand { /// Scalar type type T: RlstScalar; + #[allow(clippy::too_many_arguments)] /// Evaluate integrand for a singular quadrature rule + /// + /// # Safety + /// This method is unsafe to allow `get_unchecked` to be used unsafe fn evaluate_nonsingular( &self, test_table: &RlstArray, @@ -43,7 +47,11 @@ pub trait BoundaryIntegrand { trial_geometry: &impl CellGeometry::Real>, ) -> Self::T; + #[allow(clippy::too_many_arguments)] /// Evaluate integrand for a non-singular quadrature rule + /// + /// # Safety + /// This method is unsafe to allow `get_unchecked` to be used unsafe fn evaluate_singular( &self, test_table: &RlstArray, From 745de1737ec789bac7bcd7563e9334f37d29303f Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Wed, 21 Aug 2024 16:55:40 +0100 Subject: [PATCH 13/13] simplify --- src/assembly/boundary/assemblers.rs | 30 +++++++++-------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/src/assembly/boundary/assemblers.rs b/src/assembly/boundary/assemblers.rs index 4f2dae9b..f9ffe15d 100644 --- a/src/assembly/boundary/assemblers.rs +++ b/src/assembly/boundary/assemblers.rs @@ -151,19 +151,14 @@ where /// Assemble the contribution to the terms of a matrix for a batch of pairs of adjacent cells #[allow(clippy::too_many_arguments)] -fn assemble_batch_singular< - T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, ->( +fn assemble_batch_singular>( assembler: &impl BoundaryAssembler, deriv_size: usize, shape: [usize; 2], trial_cell_type: ReferenceCellType, test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, - test_space: &impl FunctionSpace, + trial_space: &Space, + test_space: &Space, cell_pairs: &[(usize, usize)], trial_points: &RlstArray, test_points: &RlstArray, @@ -230,20 +225,15 @@ fn assemble_batch_singular< /// Assemble the contribution to the terms of a matrix for a batch of non-adjacent cells #[allow(clippy::too_many_arguments)] -fn assemble_batch_nonadjacent< - T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, ->( +fn assemble_batch_nonadjacent>( assembler: &impl BoundaryAssembler, deriv_size: usize, output: &RawData2D, trial_cell_type: ReferenceCellType, test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, + trial_space: &Space, trial_cells: &[usize], - test_space: &impl FunctionSpace, + test_space: &Space, test_cells: &[usize], trial_points: &RlstArray, trial_weights: &[T::Real], @@ -319,17 +309,15 @@ fn assemble_batch_nonadjacent< #[allow(clippy::too_many_arguments)] fn assemble_batch_singular_correction< T: RlstScalar + MatrixInverse, - TestGrid: Grid, - TrialGrid: Grid, - Element: FiniteElement + Sync, + Space: FunctionSpace, >( assembler: &impl BoundaryAssembler, deriv_size: usize, shape: [usize; 2], trial_cell_type: ReferenceCellType, test_cell_type: ReferenceCellType, - trial_space: &impl FunctionSpace, - test_space: &impl FunctionSpace, + trial_space: &Space, + test_space: &Space, cell_pairs: &[(usize, usize)], trial_points: &RlstArray, trial_weights: &[T::Real],