Skip to content

Commit

Permalink
Refactored boundary assembler
Browse files Browse the repository at this point in the history
  • Loading branch information
@tbetcke
tbetcke committed Nov 3, 2024
1 parent f67b777 commit 65401b5
Show file tree
Hide file tree
Showing 15 changed files with 1,972 additions and 1,863 deletions.
26 changes: 19 additions & 7 deletions benches/assembly_benchmark.rs
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,11 @@
use bempp::assembly::boundary::BoundaryAssembler;
use bempp::assembly::boundary::integrands::SingleLayerBoundaryIntegrand;
use bempp::assembly::boundary::{BoundaryAssembler, BoundaryAssemblerOptions};
use bempp::assembly::kernels::KernelEvaluator;
use bempp::function::SerialFunctionSpace;
use bempp::traits::{BoundaryAssembly, FunctionSpace};
use criterion::{criterion_group, criterion_main, Criterion};
use green_kernels::laplace_3d::Laplace3dKernel;
use green_kernels::types::GreenKernelEvalType;
use ndelement::ciarlet::LagrangeElementFamily;
use ndelement::types::{Continuity, ReferenceCellType};
use ndgrid::shapes::regular_sphere;
Expand All @@ -19,21 +23,29 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
let mut matrix = rlst_dynamic_array2!(f64, [space.global_size(), space.global_size()]);

let colouring = space.cell_colouring();
let mut a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();
a.set_quadrature_degree(ReferenceCellType::Triangle, 16);
a.set_singular_quadrature_degree(

let mut assembler = BoundaryAssembler::<f64, _, _>::new(
SingleLayerBoundaryIntegrand::new(),
KernelEvaluator::new(Laplace3dKernel::new(), GreenKernelEvalType::Value),
BoundaryAssemblerOptions::default(),
1,
0,
);

assembler.set_quadrature_degree(ReferenceCellType::Triangle, 16);
assembler.set_singular_quadrature_degree(
(ReferenceCellType::Triangle, ReferenceCellType::Triangle),
4,
);
a.set_batch_size(128);
assembler.set_batch_size(128);

group.bench_function(
format!(
"Assembly of singular terms of {}x{} matrix",
space.global_size(),
space.global_size()
),
|b| b.iter(|| a.assemble_singular_into_dense(&mut matrix, &space, &space)),
|b| b.iter(|| assembler.assemble_singular_into_dense(&mut matrix, &space, &space)),
);
group.bench_function(
format!(
Expand All @@ -43,7 +55,7 @@ pub fn assembly_parts_benchmark(c: &mut Criterion) {
),
|b| {
b.iter(|| {
a.assemble_nonsingular_into_dense(
assembler.assemble_nonsingular_into_dense(
&mut matrix,
&space,
&space,
Expand Down
97 changes: 46 additions & 51 deletions examples/assembly.rs
View file
Original file line number Diff line number Diff line change
@@ -1,84 +1,79 @@
use bempp::assembly::boundary::BoundaryAssembler;
use bempp::assembly::boundary::{BoundaryAssembler, BoundaryAssemblerOptions};
use bempp::function::SerialFunctionSpace;
use bempp::laplace::assembler::laplace_single_layer;
use bempp::traits::{BoundaryAssembly, FunctionSpace};
use ndelement::ciarlet::LagrangeElementFamily;
use ndelement::types::{Continuity, ReferenceCellType};
use ndgrid::shapes::regular_sphere;
use rlst::{rlst_dynamic_array2, RandomAccessByRef};
use rlst::{rlst_dynamic_array2, RandomAccessByRef, Shape};

fn main() {
// Create a grid, family of elements, and function space
let grid = regular_sphere(0);
let element = LagrangeElementFamily::<f64>::new(1, Continuity::Standard);
let space = SerialFunctionSpace::new(&grid, &element);

// Create an array to store the assembled discrete operator
let ndofs = space.global_size();
let mut matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);

// Create an assembler for the Laplace single layer operator
let mut a = BoundaryAssembler::<f64, _, _>::new_laplace_single_layer();

// Adjust the quadrature degree for non-singular integrals on a triangle.
// This makes the integrals use a quadrature rule with 16 points
a.set_quadrature_degree(ReferenceCellType::Triangle, 16);
let mut options = BoundaryAssemblerOptions::default();
options.set_regular_quadrature_degree(ReferenceCellType::Triangle, 16);

// Adjust the quadrature degree for singular integrals on a pair ortriangles.
// This makes the integrals use a quadrature rule based on a rule on an interval with 4 points
a.set_singular_quadrature_degree(
options.set_singular_quadrature_degree(
(ReferenceCellType::Triangle, ReferenceCellType::Triangle),
4,
);

// Assemble the discrete operator
a.assemble_into_dense(&mut matrix, &space, &space);
// Assemble the single layer Laplace boundary operator.
let matrix = laplace_single_layer(&space, &space, options.clone());

// Print the entries of the matrix
println!("Lagrange single layer matrix");
for i in 0..ndofs {
for j in 0..ndofs {
for i in 0..matrix.shape()[0] {
for j in 0..matrix.shape()[1] {
print!("{:.4} ", matrix.get([i, j]).unwrap());
}
println!();
}
println!();

// Assemble just the singular part
let mut singular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
a.assemble_singular_into_dense(&mut singular_matrix, &space, &space);
println!("Lagrange single layer matrix (singular part)");
for i in 0..ndofs {
for j in 0..ndofs {
print!("{:.4} ", singular_matrix.get([i, j]).unwrap());
}
println!();
}
println!();
// // Assemble just the singular part
// let mut singular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
// a.assemble_singular_into_dense(&mut singular_matrix, &space, &space);
// println!("Lagrange single layer matrix (singular part)");
// for i in 0..ndofs {
// for j in 0..ndofs {
// print!("{:.4} ", singular_matrix.get([i, j]).unwrap());
// }
// println!();
// }
// println!();

// For grids with a larger number of cells, the singular part will be sparse. It can be assembled as a CSR matrix as follows
println!("Lagrange single layer matrix (singular part) as CSR matrix");
let singular_sparse_matrix = a.assemble_singular_into_csr(&space, &space);
println!("indices: {:?}", singular_sparse_matrix.indices());
println!("indptr: {:?}", singular_sparse_matrix.indptr());
println!("data: {:?}", singular_sparse_matrix.data());
println!();
// // For grids with a larger number of cells, the singular part will be sparse. It can be assembled as a CSR matrix as follows
// println!("Lagrange single layer matrix (singular part) as CSR matrix");
// let singular_sparse_matrix = a.assemble_singular_into_csr(&space, &space);
// println!("indices: {:?}", singular_sparse_matrix.indices());
// println!("indptr: {:?}", singular_sparse_matrix.indptr());
// println!("data: {:?}", singular_sparse_matrix.data());
// println!();

// Assemble just the non-singular part
let colouring = space.cell_colouring();
let mut nonsingular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
a.assemble_nonsingular_into_dense(
&mut nonsingular_matrix,
&space,
&space,
&colouring,
&colouring,
);
println!("Lagrange single layer matrix (nonsingular part)");
for i in 0..ndofs {
for j in 0..ndofs {
print!("{:.4} ", nonsingular_matrix.get([i, j]).unwrap());
}
println!();
}
println!();
// // Assemble just the non-singular part
// let colouring = space.cell_colouring();
// let mut nonsingular_matrix = rlst_dynamic_array2!(f64, [ndofs, ndofs]);
// a.assemble_nonsingular_into_dense(
// &mut nonsingular_matrix,
// &space,
// &space,
// &colouring,
// &colouring,
// );
// println!("Lagrange single layer matrix (nonsingular part)");
// for i in 0..ndofs {
// for j in 0..ndofs {
// print!("{:.4} ", nonsingular_matrix.get([i, j]).unwrap());
// }
// println!();
// }
// println!();
}
Loading

0 comments on commit 65401b5

Please sign in to comment.