Skip to content

Commit

Permalink
Lots of small fixes
Browse files Browse the repository at this point in the history
  • Loading branch information
@tbetcke
tbetcke committed Jan 26, 2025
1 parent 0a30d4a commit 9f970a5
Show file tree
Hide file tree
Showing 6 changed files with 193 additions and 420 deletions.
12 changes: 12 additions & 0 deletions examples/function_evaluator.rs
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,12 @@
//! Demo the evaluation of functions in Bempp-rs
fn main() {
let universe = mpi::initialize().unwrap();
let world = universe.world();

// We create a sphere grid

let grid = bempp::shapes::regular_sphere::<f64, _>(5, 1, &world);

// We have the grid.
}
77 changes: 15 additions & 62 deletions examples/laplace_evaluator.rs
View file
Original file line number Diff line number Diff line change
@@ -1,7 +1,5 @@
//! This file implements an example Laplace evaluator test the different involved operators.
use std::rc::Rc;

use bempp::{
boundary_assemblers::BoundaryAssemblerOptions,
evaluator_tools::NeighbourEvaluator,
Expand All @@ -10,14 +8,10 @@ use bempp::{
use green_kernels::laplace_3d::Laplace3dKernel;
use mpi::traits::Communicator;
use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
use ndgrid::{
traits::{Entity, GeometryMap, Grid, ParallelGrid},
types::Ownership,
};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use rlst::{
operator::{interface::DistributedArrayVectorSpace, zero_element, Operator},
operator::{zero_element, Operator},
rlst_dynamic_array1, AsApply, MultInto, OperatorBase,
};

Expand All @@ -32,10 +26,6 @@ fn main() {
let quad_degree = 6;
// Get the number of cells in the grid.

let n_cells = grid.local_grid().entity_iter(2).count();

println!("Number of cells: {}", n_cells);

let space = bempp::function::FunctionSpace::new(
&grid,
&LagrangeElementFamily::<f64>::new(1, ndelement::types::Continuity::Discontinuous),
Expand All @@ -54,25 +44,8 @@ fn main() {

// Now let's build an evaluator.

//First initialise the index layouts.

let space_layout = Rc::new(bempp_distributed_tools::IndexLayout::from_local_counts(
space.local_space().global_size(),
&world,
));

let point_layout = Rc::new(bempp_distributed_tools::IndexLayout::from_local_counts(
quad_degree * n_cells,
&world,
));

// Instantiate function spaces.

let array_function_space =
DistributedArrayVectorSpace::<_, f64>::from_index_layout(space_layout.clone());
let point_function_space =
DistributedArrayVectorSpace::<_, f64>::from_index_layout(point_layout.clone());

let qrule = bempp_quadrature::simplex_rules::simplex_rule_triangle(quad_degree).unwrap();

let space_to_point =
Expand All @@ -83,46 +56,26 @@ fn main() {

// We now have to get all the points from the grid. We do this by iterating through all cells.

let mut points = vec![0 as f64; 3 * quad_degree * n_cells];

let geometry_map = grid
.local_grid()
.geometry_map(ReferenceCellType::Triangle, &qrule.points);

for cell in grid
.local_grid()
.entity_iter(2)
.filter(|e| matches!(e.ownership(), Ownership::Owned))
{
let start_index = 3 * point_function_space
.index_layout()
.global2local(world.rank() as usize, qrule.npoints * cell.global_index())
.unwrap();
geometry_map.points(
cell.local_index(),
&mut points[start_index..start_index + 3 * qrule.npoints],
);
}

// let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
// &points,
// &points,
// green_kernels::types::GreenKernelEvalType::Value,
// true,
// Laplace3dKernel::default(),
// &world,
// );
let points = bempp::evaluator_tools::grid_points_from_space(&space, &qrule.points);

let kernel_evaluator = bempp::greens_function_evaluators::kifmm_evaluator::KiFmmEvaluator::new(
&points, &points, 1, 3, 5, &world,
let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
&points,
&points,
green_kernels::types::GreenKernelEvalType::Value,
true,
Laplace3dKernel::default(),
&world,
);

// let kernel_evaluator = bempp::greens_function_evaluators::kifmm_evaluator::KiFmmEvaluator::new(
// &points, &points, 1, 3, 5, &world,
// );

let correction = NeighbourEvaluator::new(
&space,
&qrule.points,
Laplace3dKernel::default(),
green_kernels::types::GreenKernelEvalType::Value,
space.local_space().support_cells(),
&grid,
);

let corrected_evaluator = kernel_evaluator.r().sum(correction.r().scale(-1.0));
Expand All @@ -133,7 +86,7 @@ fn main() {

let laplace_evaluator = (point_to_space.product(prod1)).sum(singular_operator.r());

let mut x = zero_element(array_function_space.clone());
let mut x = zero_element(laplace_evaluator.domain());
x.view_mut()
.local_mut()
.fill_from_equally_distributed(&mut rng);
Expand Down
77 changes: 23 additions & 54 deletions examples/neighbour_evaluator.rs
View file
Original file line number Diff line number Diff line change
@@ -1,79 +1,48 @@
//! Test the neighbour evaluator
use std::rc::Rc;

use bempp::evaluator_tools::NeighbourEvaluator;
use green_kernels::laplace_3d::Laplace3dKernel;
use itertools::Itertools;
use mpi::traits::Communicator;
use ndelement::types::ReferenceCellType;
use ndgrid::{
traits::{Entity, GeometryMap, Grid, ParallelGrid},
types::Ownership,
};
use rand::SeedableRng;
use rand_chacha::ChaCha8Rng;
use rlst::{
operator::{interface::DistributedArrayVectorSpace, zero_element},
rlst_dynamic_array2, AsApply, IndexLayout, RandomAccessMut, RawAccess,
use bempp::{
boundary_assemblers::BoundaryAssemblerOptions,
evaluator_tools::{grid_points_from_space, NeighbourEvaluator},
};
use green_kernels::laplace_3d::Laplace3dKernel;
use ndelement::{ciarlet::LagrangeElementFamily, types::ReferenceCellType};
use rlst::{operator::zero_element, AsApply, OperatorBase, RandomAccessMut};

fn main() {
let universe = mpi::initialize().unwrap();
let world = universe.world();

let n_points = 5;

let mut rng = ChaCha8Rng::seed_from_u64(world.rank() as u64);
let grid = bempp::shapes::regular_sphere::<f64, _>(5, 1, &world);

let mut points = rlst_dynamic_array2!(f64, [2, n_points]);
points.fill_from_equally_distributed(&mut rng);
let quad_degree = 6;
// Get the number of cells in the grid.

let grid = bempp::shapes::regular_sphere::<f64, _>(3, 1, &world);

// Now get the active cells on the current process.
let space = bempp::function::FunctionSpace::new(
&grid,
&LagrangeElementFamily::<f64>::new(1, ndelement::types::Continuity::Discontinuous),
);

let n_cells = grid
.local_grid()
.entity_iter(2)
.filter(|e| matches!(e.ownership(), Ownership::Owned))
.count();
let mut options = BoundaryAssemblerOptions::default();
options.set_regular_quadrature_degree(ReferenceCellType::Triangle, quad_degree);

let index_layout = Rc::new(IndexLayout::from_local_counts(n_cells * n_points, &world));
let quad_degree = options
.get_regular_quadrature_degree(ReferenceCellType::Triangle)
.unwrap();

let space = DistributedArrayVectorSpace::<_, f64>::from_index_layout(index_layout.clone());
let qrule = bempp_quadrature::simplex_rules::simplex_rule_triangle(quad_degree).unwrap();

let neighbour_evaluator = NeighbourEvaluator::new(
points.data(),
&space,
&qrule.points,
Laplace3dKernel::default(),
green_kernels::types::GreenKernelEvalType::Value,
&(0..grid.local_grid().entity_count(ReferenceCellType::Triangle)).collect_vec(),
&grid,
);

// We now manually test the evaluator. For that we first create a dense evaluator so that we have comparison.

// For the evaluator we need all the points.

let mut physical_points = vec![0 as f64; 3 * n_points * n_cells];

let geometry_map = grid
.local_grid()
.geometry_map(ReferenceCellType::Triangle, points.data());

for cell in grid
.local_grid()
.entity_iter(2)
.filter(|e| matches!(e.ownership(), Ownership::Owned))
{
let start_index = 3 * index_layout
.global2local(world.rank() as usize, n_points * cell.global_index())
.unwrap();
geometry_map.points(
cell.local_index(),
&mut physical_points[start_index..start_index + 3 * n_points],
);
}
let physical_points = grid_points_from_space(&space, &qrule.points);

let kernel_evaluator = bempp::greens_function_evaluators::dense_evaluator::DenseEvaluator::new(
&physical_points,
Expand All @@ -84,7 +53,7 @@ fn main() {
&world,
);

let mut x = zero_element(space.clone());
let mut x = zero_element(kernel_evaluator.domain());

*x.view_mut().local_mut().get_mut([0]).unwrap() = 1.0;
*x.view_mut().local_mut().get_mut([1]).unwrap() = 2.0;
Expand Down
Loading

0 comments on commit 9f970a5

Please sign in to comment.