Skip to content

Commit

Permalink
wrap boundary assembler options
Browse files Browse the repository at this point in the history
  • Loading branch information
mscroggs committed Sep 18, 2024
1 parent d57f834 commit 0e32fa3
Show file tree
Hide file tree
Showing 15 changed files with 1,894 additions and 15 deletions.
4 changes: 2 additions & 2 deletions examples/assembly.rs
Original file line number Diff line number Diff line change
Expand Up @@ -21,11 +21,11 @@ fn main() {

// Adjust the quadrature degree for non-singular integrals on a triangle.
// This makes the integrals use a quadrature rule with 16 points
a.quadrature_degree(ReferenceCellType::Triangle, 16);
a.set_quadrature_degree(ReferenceCellType::Triangle, 16);

// Adjust the quadrature degree for singular integrals on a pair ortriangles.
// This makes the integrals use a quadrature rule based on a rule on an interval with 4 points
a.singular_quadrature_degree(
a.set_singular_quadrature_degree(
(ReferenceCellType::Triangle, ReferenceCellType::Triangle),
4,
);
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,7 @@ dependencies = [
"ndelement",
"ndgrid"
]
packages = ["bempp"]
packages = ["bempp", "bempp.assembly"]

[project.urls]
homepage = "https://github.com/bempp/bempp-rs"
Expand Down
2 changes: 1 addition & 1 deletion python/bempp/__init__.py
Original file line number Diff line number Diff line change
@@ -1,3 +1,3 @@
"""Bempp."""

from bempp import function_space
from bempp import function_space, assembly
3 changes: 3 additions & 0 deletions python/bempp/assembly/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,3 @@
"""Assembly."""

from bempp.assembly import boundary
112 changes: 112 additions & 0 deletions python/bempp/assembly/boundary.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,112 @@
"""Boundary assembly."""

import typing
import numpy as np
from bempp._bempprs import lib as _lib
from ndelement.reference_cell import ReferenceCellType
from enum import Enum
from _cffi_backend import _CDataBase

_dtypes = {
0: np.float32,
1: np.float64,
}
_ctypes = {
np.float32: "float",
np.float64: "double",
}


def _dtype_value(dtype):
"""Get the u8 identifier for a data type."""
for i, j in _dtypes.items():
if j == dtype:
return i
raise TypeError("Invalid data type")


class OperatorType(Enum):
"""Operator type."""

SingleLayer = 0
DoubleLayer = 1
AdjointDoubleLayer = 2
Hypersingular = 3
ElectricField = 4
MagneticField = 5


class BoundaryAssembler(object):
"""Boundary assembler."""

def __init__(self, rs_assembler: _CDataBase, owned: bool = True):
"""Initialise."""
self._rs_assembler = rs_assembler
self._owned = owned

def __del__(self):
"""Delete."""
if self._owned:
_lib.free_boundary_assembler(self._rs_assembler)

def set_quadrature_degree(self, cell: ReferenceCellType, degree: int):
"""Set the (non-singular) quadrature degree for a cell type."""
if not _lib.boundary_assembler_has_quadrature_degree(self._rs_assembler, cell.value):
raise ValueError(f"Invalid cell type: {cell}")
_lib.boundary_assembler_set_quadrature_degree(self._rs_assembler, cell.value, degree)

def quadrature_degree(self, cell: ReferenceCellType) -> int:
"""Get the (non-singular) quadrature degree for a cell type."""
if not _lib.boundary_assembler_has_quadrature_degree(self._rs_assembler, cell.value):
raise ValueError(f"Invalid cell type: {cell}")
return _lib.boundary_assembler_quadrature_degree(self._rs_assembler, cell.value)

def set_singular_quadrature_degree(
self, cell0: ReferenceCellType, cell1: ReferenceCellType, degree: int
):
"""Set the singular quadrature degree for a cell type."""
if not _lib.boundary_assembler_has_singular_quadrature_degree(
self._rs_assembler, cell0.value, cell1.value
):
raise ValueError(f"Invalid cell pair: {cell0} and {cell1}")
_lib.boundary_assembler_set_singular_quadrature_degree(
self._rs_assembler, cell0.value, cell1.value, degree
)

def singular_quadrature_degree(self, cell0: ReferenceCellType, cell1: ReferenceCellType) -> int:
"""Get the singular_quadrature degree for a cell type."""
if not _lib.boundary_assembler_has_singular_quadrature_degree(
self._rs_assembler, cell0.value, cell1.value
):
raise ValueError(f"Invalid cell pair: {cell0} and {cell1}")
return _lib.boundary_assembler_singular_quadrature_degree(
self._rs_assembler, cell0.value, cell1.value
)

def set_batch_size(self, batch_size: int):
"""Set the batch size."""
_lib.boundary_assembler_set_batch_size(self._rs_assembler, batch_size)

@property
def batch_size(self) -> int:
"""Get the batch size."""
return _lib.boundary_assembler_batch_size(self._rs_assembler)

@property
def dtype(self):
"""Data type."""
return _dtypes[_lib.boundary_assembler_dtype(self._rs_assembler)]

@property
def _ctype(self):
"""C data type."""
return _ctypes[self.dtype]


def create_laplace_assembler(
operator_type: OperatorType, dtype: typing.Type[np.floating] = np.float64
):
"""Create a Laplace assembler."""
return BoundaryAssembler(
_lib.laplace_boundary_assembler_new(operator_type.value, _dtype_value(dtype))
)
7 changes: 5 additions & 2 deletions python/bempp/function_space.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
from ndgrid.ownership import Owned, Ghost
from ndelement.ciarlet import ElementFamily, CiarletElement
from ndelement.reference_cell import ReferenceCellType
from _cffi_backend import _CDataBase

_dtypes = {
0: np.float32,
Expand All @@ -23,13 +24,15 @@
class FunctionSpace(object):
"""Function space."""

def __init__(self, rs_space):
def __init__(self, rs_space: _CDataBase, owned: bool = True):
"""Initialise."""
self._rs_space = rs_space
self._owned = True

def __del__(self):
"""Delete."""
_lib.free_space(self._rs_space)
if self._owned:
_lib.free_space(self._rs_space)

def element(self, entity: ReferenceCellType) -> CiarletElement:
"""Get the grid that this space is defined on."""
Expand Down
37 changes: 37 additions & 0 deletions python/test/test_boundary_assembler.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
import pytest
import numpy as np
from ndelement.reference_cell import ReferenceCellType
from bempp.assembly.boundary import OperatorType, create_laplace_assembler


@pytest.mark.parametrize(
"otype",
[
OperatorType.SingleLayer,
OperatorType.DoubleLayer,
OperatorType.AdjointDoubleLayer,
OperatorType.Hypersingular,
],
)
def test_create_assembler(otype):
a = create_laplace_assembler(otype)

assert a.quadrature_degree(ReferenceCellType.Triangle) != 3
a.set_quadrature_degree(ReferenceCellType.Triangle, 3)
assert a.quadrature_degree(ReferenceCellType.Triangle) == 3
with pytest.raises(ValueError):
a.set_quadrature_degree(ReferenceCellType.Interval, 3)
with pytest.raises(ValueError):
a.quadrature_degree(ReferenceCellType.Interval)

assert a.singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle) != 3
a.set_singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle, 3)
assert a.singular_quadrature_degree(ReferenceCellType.Triangle, ReferenceCellType.Triangle) == 3
with pytest.raises(ValueError):
a.set_singular_quadrature_degree(ReferenceCellType.Interval, ReferenceCellType.Interval, 3)

assert a.batch_size != 4
a.set_batch_size(4)
assert a.batch_size == 4

assert a.dtype == np.float64
2 changes: 1 addition & 1 deletion src/assembly/boundary.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
//! Assembly of boundary operators
mod assemblers;
pub(crate) mod cell_pair_assemblers;
mod integrands;
pub mod integrands;

pub use assemblers::BoundaryAssembler;

Expand Down
27 changes: 24 additions & 3 deletions src/assembly/boundary/assemblers.rs
Original file line number Diff line number Diff line change
Expand Up @@ -484,12 +484,17 @@ impl<
}

/// Set (non-singular) quadrature degree for a cell type
pub fn quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
pub fn set_quadrature_degree(&mut self, cell: ReferenceCellType, degree: usize) {
*self.options.quadrature_degrees.get_mut(&cell).unwrap() = degree;
}

/// Get (non-singular) quadrature degree for a cell type
pub fn quadrature_degree(&self, cell: ReferenceCellType) -> Option<usize> {
self.options.quadrature_degrees.get(&cell).map(|i| *i)

Check failure on line 493 in src/assembly/boundary/assemblers.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,strict")

you are using an explicit closure for copying elements

Check failure on line 493 in src/assembly/boundary/assemblers.rs

View workflow job for this annotation

GitHub Actions / Rust style checks

you are using an explicit closure for copying elements
}

/// Set singular quadrature degree for a pair of cell types
pub fn singular_quadrature_degree(
pub fn set_singular_quadrature_degree(
&mut self,
cells: (ReferenceCellType, ReferenceCellType),
degree: usize,
Expand All @@ -501,11 +506,27 @@ impl<
.unwrap() = degree;
}

/// Get singular quadrature degree for a pair of cell types
pub fn singular_quadrature_degree(
&self,
cells: (ReferenceCellType, ReferenceCellType),
) -> Option<usize> {
self.options

Check failure on line 514 in src/assembly/boundary/assemblers.rs

View workflow job for this annotation

GitHub Actions / Rust style checks (--features "mpi,strict")

you are using an explicit closure for copying elements

Check failure on line 514 in src/assembly/boundary/assemblers.rs

View workflow job for this annotation

GitHub Actions / Rust style checks

you are using an explicit closure for copying elements
.singular_quadrature_degrees
.get(&cells)
.map(|i| *i)
}

/// Set the maximum size of a batch of cells to send to an assembly function
pub fn batch_size(&mut self, size: usize) {
pub fn set_batch_size(&mut self, size: usize) {
self.options.batch_size = size;
}

/// Get the maximum size of a batch of cells to send to an assembly function
pub fn batch_size(&self) -> usize {
self.options.batch_size
}

/// Assemble the singular contributions
fn assemble_singular<Space: FunctionSpace<T = T> + Sync>(
&self,
Expand Down
2 changes: 2 additions & 0 deletions src/assembly/boundary/integrands.rs
Original file line number Diff line number Diff line change
Expand Up @@ -27,6 +27,7 @@ pub struct BoundaryIntegrandSum<
impl<T: RlstScalar, I0: BoundaryIntegrand<T = T>, I1: BoundaryIntegrand<T = T>>
BoundaryIntegrandSum<T, I0, I1>
{
/// Create new
pub fn new(integrand0: I0, integrand1: I1) -> Self {
Self {
integrand0,
Expand Down Expand Up @@ -63,6 +64,7 @@ pub struct BoundaryIntegrandScalarProduct<T: RlstScalar, I: BoundaryIntegrand<T
}

impl<T: RlstScalar, I: BoundaryIntegrand<T = T>> BoundaryIntegrandScalarProduct<T, I> {
/// Create new
pub fn new(scalar: T, integrand: I) -> Self {
Self { scalar, integrand }
}
Expand Down
2 changes: 2 additions & 0 deletions src/assembly/boundary/integrands/adjoint_double_layer.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
use rlst::RlstScalar;

/// Integrand for an adjoint double layer boundary operator
pub struct AdjointDoubleLayerBoundaryIntegrand<T: RlstScalar> {
_t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> AdjointDoubleLayerBoundaryIntegrand<T> {
/// Create new
pub fn new() -> Self {
Self {
_t: std::marker::PhantomData,
Expand Down
2 changes: 2 additions & 0 deletions src/assembly/boundary/integrands/double_layer.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
use rlst::RlstScalar;

/// Integrand for a double layer boundary operator
pub struct DoubleLayerBoundaryIntegrand<T: RlstScalar> {
_t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> DoubleLayerBoundaryIntegrand<T> {
/// Create new
pub fn new() -> Self {
Self {
_t: std::marker::PhantomData,
Expand Down
4 changes: 4 additions & 0 deletions src/assembly/boundary/integrands/hypersingular.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
use rlst::RlstScalar;

/// Integrand for the curl curl term of a hypersingular boundary operator
pub struct HypersingularCurlCurlBoundaryIntegrand<T: RlstScalar> {
_t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> HypersingularCurlCurlBoundaryIntegrand<T> {
/// Create new
pub fn new() -> Self {
Self {
_t: std::marker::PhantomData,
Expand Down Expand Up @@ -52,11 +54,13 @@ unsafe impl<T: RlstScalar> BoundaryIntegrand for HypersingularCurlCurlBoundaryIn
}
}

/// Integrand for the normal normal term of a hypersingular boundary operator
pub struct HypersingularNormalNormalBoundaryIntegrand<T: RlstScalar> {
_t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> HypersingularNormalNormalBoundaryIntegrand<T> {
/// Create new
pub fn new() -> Self {
Self {
_t: std::marker::PhantomData,
Expand Down
2 changes: 2 additions & 0 deletions src/assembly/boundary/integrands/single_layer.rs
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,13 @@
use crate::traits::{Access1D, Access2D, BoundaryIntegrand, GeometryAccess};
use rlst::RlstScalar;

/// Integrand for a single layer boundary operator
pub struct SingleLayerBoundaryIntegrand<T: RlstScalar> {
_t: std::marker::PhantomData<T>,
}

impl<T: RlstScalar> SingleLayerBoundaryIntegrand<T> {
/// Create new
pub fn new() -> Self {
Self {
_t: std::marker::PhantomData,
Expand Down
Loading

0 comments on commit 0e32fa3

Please sign in to comment.