Skip to content

Commit a6cd0cb

Browse files
w1th0utnam3w1th0utnam3
committed
Reduce code duplication
1 parent d14f817 commit a6cd0cb

File tree

1 file changed

+64
-129
lines changed

1 file changed

+64
-129
lines changed

‎splashsurf_lib/src/dense_subdomains.rs‎

Lines changed: 64 additions & 129 deletions
Original file line numberDiff line numberDiff line change
@@ -657,6 +657,62 @@ pub(crate) struct SurfacePatch<I: Index, R: Real> {
657657
pub exterior_vertex_edge_indices: Vec<(I, EdgeIndex<I>)>,
658658
}
659659

660+
/// Returns the lower and upper subdomain grid vertex indices possibly influenced by the given particle (upper index is exclusive)
661+
#[inline(always)]
662+
fn particle_influence_aabb<I: Index, R: Real>(
663+
particle: &Vector3<R>,
664+
subdomain_mc_grid: &UniformCartesianCubeGrid3d<I, R>,
665+
cube_radius: I,
666+
) -> ([I; 3], [I; 3]) {
667+
let extents = subdomain_mc_grid.points_per_dim();
668+
669+
// Note: This assumes that enclosing_cell can return negative indices for ghost particles
670+
// (they are outside the subdomain grid but part of the particle list).
671+
672+
// Get grid cell containing particle
673+
let particle_cell = subdomain_mc_grid.enclosing_cell(particle);
674+
675+
// Compute lower and upper bounds of the grid points possibly affected by the particle
676+
// We want to loop over the vertices of the enclosing cells plus all points in `cube_radius` distance from the cell
677+
678+
let lower = [0, 1, 2].map(|i| {
679+
particle_cell[i]
680+
.saturating_sub(&cube_radius)
681+
.max(I::zero())
682+
.min(extents[i])
683+
});
684+
685+
let upper = [0, 1, 2].map(|i| {
686+
// We add 2 because
687+
// - we want to loop over all grid points of the cell (+1 for upper points) + the radius
688+
// - the upper range limit is exclusive (+1)
689+
(particle_cell[i] + cube_radius + I::two())
690+
.min(extents[i])
691+
.max(I::zero())
692+
});
693+
694+
(lower, upper)
695+
}
696+
697+
/// Converts a local (subdomain) grid point index to a global grid point index
698+
fn local_to_global_point_ijk<I: Index>(
699+
local_point_ijk: [I; 3],
700+
subdomain_ijk: [I; 3],
701+
cells_per_subdomain: [I; 3],
702+
) -> [GlobalIndex; 3] {
703+
let local_point_ijk = local_point_ijk.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
704+
let subdomain_ijk = subdomain_ijk.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
705+
let cells_per_subdomain =
706+
cells_per_subdomain.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
707+
let [i, j, k] = local_point_ijk;
708+
709+
[
710+
subdomain_ijk[0] * cells_per_subdomain[0] + i,
711+
subdomain_ijk[1] * cells_per_subdomain[1] + j,
712+
subdomain_ijk[2] * cells_per_subdomain[2] + k,
713+
]
714+
}
715+
660716
/// Auto-dispatching density grid loop for f32/i64: chooses AVX2+FMA on x86, NEON on aarch64, otherwise scalar fallback
661717
pub fn density_grid_loop_auto<K: SymmetricKernel3d<f32>>(
662718
levelset_grid: &mut [f32],
@@ -741,36 +797,13 @@ pub fn density_grid_loop_scalar<I: Index, R: Real, K: SymmetricKernel3d<R>>(
741797
) {
742798
profile!("density grid loop (scalar)");
743799
let mc_grid = subdomain_mc_grid;
744-
let extents = mc_grid.points_per_dim();
745800

746801
for (p_i, rho_i) in subdomain_particles
747802
.iter()
748803
.copied()
749804
.zip(subdomain_particle_densities.iter().copied())
750805
{
751-
// Note: this loop assumes that enclosing_cell can return negative indices for ghost particles
752-
753-
// Get grid cell containing particle
754-
let particle_cell = mc_grid.enclosing_cell(&p_i);
755-
756-
// Compute lower and upper bounds of the grid points possibly affected by the particle
757-
// We want to loop over the vertices of the enclosing cells plus all points in `cube_radius` distance from the cell
758-
759-
let lower = [0, 1, 2].map(|i| {
760-
particle_cell[i]
761-
.saturating_sub(&cube_radius)
762-
.max(I::zero())
763-
.min(extents[i])
764-
});
765-
766-
let upper = [0, 1, 2].map(|i| {
767-
// We add 2 because
768-
// - we want to loop over all grid points of the cell (+1 for upper points) + the radius
769-
// - the upper range limit is exclusive (+1)
770-
(particle_cell[i] + cube_radius + I::two())
771-
.min(extents[i])
772-
.max(I::zero())
773-
});
806+
let (lower, upper) = particle_influence_aabb(&p_i, &mc_grid, cube_radius);
774807

775808
// Loop over all grid points around the enclosing cell
776809
for i in I::range(lower[0], upper[0]).iter() {
@@ -783,26 +816,6 @@ pub fn density_grid_loop_scalar<I: Index, R: Real, K: SymmetricKernel3d<R>>(
783816

784817
let mc_cells_per_subdomain = mc_grid.cells_per_dim();
785818

786-
fn local_to_global_point_ijk<I: Index>(
787-
local_point_ijk: [I; 3],
788-
subdomain_ijk: [I; 3],
789-
cells_per_subdomain: [I; 3],
790-
) -> [GlobalIndex; 3] {
791-
let local_point_ijk =
792-
local_point_ijk.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
793-
let subdomain_ijk =
794-
subdomain_ijk.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
795-
let cells_per_subdomain =
796-
cells_per_subdomain.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
797-
let [i, j, k] = local_point_ijk;
798-
799-
[
800-
subdomain_ijk[0] * cells_per_subdomain[0] + i,
801-
subdomain_ijk[1] * cells_per_subdomain[1] + j,
802-
subdomain_ijk[2] * cells_per_subdomain[2] + k,
803-
]
804-
}
805-
806819
// Use global coordinate calculation for consistency with neighboring domains
807820
let global_point_ijk = local_to_global_point_ijk(
808821
point_ijk,
@@ -888,29 +901,9 @@ pub fn density_grid_loop_neon<K: SymmetricKernel3d<f32>>(
888901
{
889902
let v_i = particle_rest_mass / rho_i;
890903

891-
// Note: this loop assumes that enclosing_cell can return negative indices for ghost particles
892-
893-
// Get grid cell containing particle
894-
let particle_cell = mc_grid.enclosing_cell(&p_i);
895-
896-
// Compute lower and upper bounds of the grid points possibly affected by the particle
897-
// We want to loop over the vertices of the enclosing cells plus all points in `cube_radius` distance from the cell
898-
899-
let lower = [0, 1, 2].map(|i| {
900-
particle_cell[i]
901-
.saturating_sub(cube_radius)
902-
.max(0)
903-
.min(mc_points[i]) as usize
904-
});
905-
906-
let upper = [0, 1, 2].map(|i| {
907-
// We add 2 because
908-
// - we want to loop over all grid points of the cell (+1 for upper points) + the radius
909-
// - the upper range limit is exclusive (+1)
910-
(particle_cell[i] + cube_radius + 2)
911-
.min(mc_points[i])
912-
.max(0) as usize
913-
});
904+
let (lower, upper) = particle_influence_aabb(&p_i, &mc_grid, cube_radius);
905+
let lower = lower.map(|i| i as usize);
906+
let upper = upper.map(|i| i as usize);
914907

915908
let remainder = (upper[2] - lower[2]) % LANES;
916909
let upper_k_aligned = upper[2] - remainder;
@@ -1052,21 +1045,9 @@ pub fn density_grid_loop_avx<K: SymmetricKernel3d<f32>>(
10521045
{
10531046
let v_i = particle_rest_mass / rho_i;
10541047

1055-
// Get grid cell containing particle
1056-
let particle_cell = mc_grid.enclosing_cell(&p_i);
1057-
1058-
// Bounds of possibly affected grid points
1059-
let lower = [0, 1, 2].map(|ii| {
1060-
particle_cell[ii]
1061-
.saturating_sub(cube_radius)
1062-
.max(0)
1063-
.min(mc_points[ii]) as usize
1064-
});
1065-
let upper = [0, 1, 2].map(|ii| {
1066-
(particle_cell[ii] + cube_radius + 2)
1067-
.min(mc_points[ii])
1068-
.max(0) as usize
1069-
});
1048+
let (lower, upper) = particle_influence_aabb(&p_i, &mc_grid, cube_radius);
1049+
let lower = lower.map(|i| i as usize);
1050+
let upper = upper.map(|i| i as usize);
10701051

10711052
let remainder = (upper[2] - lower[2]) % LANES;
10721053
let upper_k_aligned = upper[2] - remainder;
@@ -1531,37 +1512,12 @@ pub(crate) fn reconstruction<I: Index, R: Real>(
15311512

15321513
{
15331514
profile!("density grid loop");
1534-
1535-
let extents = mc_grid.points_per_dim();
1536-
15371515
for (p_i, rho_i) in subdomain_particles
15381516
.iter()
15391517
.copied()
15401518
.zip(subdomain_particle_densities.iter().copied())
15411519
{
1542-
// Note: this loop assumes that enclosing_cell can return negative indices for ghost particles
1543-
1544-
// Get grid cell containing particle
1545-
let particle_cell = mc_grid.enclosing_cell(&p_i);
1546-
1547-
// Compute lower and upper bounds of the grid points possibly affected by the particle
1548-
// We want to loop over the vertices of the enclosing cells plus all points in `cube_radius` distance from the cell
1549-
1550-
let lower = [0, 1, 2].map(|i| {
1551-
particle_cell[i]
1552-
.saturating_sub(&cube_radius)
1553-
.max(I::zero())
1554-
.min(extents[i])
1555-
});
1556-
1557-
let upper = [0, 1, 2].map(|i| {
1558-
// We add 2 because
1559-
// - we want to loop over all grid points of the cell (+1 for upper points) + the radius
1560-
// - the upper range limit is exclusive (+1)
1561-
(particle_cell[i] + cube_radius + I::two())
1562-
.min(extents[i])
1563-
.max(I::zero())
1564-
});
1520+
let (lower, upper) = particle_influence_aabb(&p_i, &mc_grid, cube_radius);
15651521

15661522
// Loop over all grid points around the enclosing cell
15671523
for i in I::range(lower[0], upper[0]).iter() {
@@ -1571,31 +1527,10 @@ pub(crate) fn reconstruction<I: Index, R: Real>(
15711527
let local_point = mc_grid
15721528
.get_point(point_ijk)
15731529
.expect("point has to be part of the subdomain grid");
1574-
//let point_coordinates = mc_grid.point_coordinates(&point);
15751530

15761531
let subdomain_ijk = subdomain_idx.index();
15771532
let mc_cells_per_subdomain = mc_grid.cells_per_dim();
15781533

1579-
fn local_to_global_point_ijk<I: Index>(
1580-
local_point_ijk: [I; 3],
1581-
subdomain_ijk: [I; 3],
1582-
cells_per_subdomain: [I; 3],
1583-
) -> [GlobalIndex; 3] {
1584-
let local_point_ijk = local_point_ijk
1585-
.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
1586-
let subdomain_ijk = subdomain_ijk
1587-
.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
1588-
let cells_per_subdomain = cells_per_subdomain
1589-
.map(|i| <GlobalIndex as NumCast>::from(i).unwrap());
1590-
let [i, j, k] = local_point_ijk;
1591-
1592-
[
1593-
subdomain_ijk[0] * cells_per_subdomain[0] + i,
1594-
subdomain_ijk[1] * cells_per_subdomain[1] + j,
1595-
subdomain_ijk[2] * cells_per_subdomain[2] + k,
1596-
]
1597-
}
1598-
15991534
// Use global coordinate calculation for consistency with neighboring domains
16001535
let global_point_ijk = local_to_global_point_ijk(
16011536
point_ijk,

0 commit comments

Comments
 (0)