lindhard/geometry/
voxel.rs

1//! The 3D [`VoxelGrid`] with its exact DDA ray traversal. See the
2//! [module docs](super) for the origin, flattening order, face ownership,
3//! boundary behaviour and the traversal algorithm.
4
5use super::{Exit, ExitOutcome, Face, Flight, Geometry, GeometryError};
6use crate::material::Material;
7
8/// What a voxel grid does to a particle that reaches the ends of one axis.
9#[derive(Debug, Clone, Copy, PartialEq, Eq)]
10pub enum Boundary {
11    /// The particle leaves the target (vacuum beyond the face).
12    Vacuum,
13    /// The particle re-enters at the opposite face, unchanged in direction
14    /// and energy.
15    Periodic,
16}
17
18/// A regular 3D grid of voxels, each filled by one material. See the module
19/// docs for the origin, flattening order, face ownership, boundary
20/// behaviour and the traversal algorithm.
21///
22/// ```
23/// use lindhard::geometry::{Boundary, Geometry, VoxelGrid};
24/// use lindhard::material::Material;
25///
26/// let si = Material::from_atom_fractions(&[(14, 1.0)], None).unwrap();
27/// let c = Material::from_atom_fractions(&[(6, 1.0)], None).unwrap();
28/// // 2 x 1 x 1 voxels of 10 nm: silicon in front of carbon.
29/// let grid = VoxelGrid::new(
30///     vec![si, c],
31///     [2, 1, 1],
32///     [1e-8; 3],
33///     vec![0, 1],
34///     [Boundary::Vacuum, Boundary::Periodic, Boundary::Periodic],
35/// )
36/// .unwrap();
37/// assert_eq!(grid.locate([1.5e-8, 0.5e-8, 0.5e-8]), Some(1));
38/// ```
39#[derive(Debug, Clone, PartialEq)]
40pub struct VoxelGrid {
41    materials: Vec<Material>,
42    dims: [usize; 3],
43    spacing_m: [f64; 3],
44    extent_m: [f64; 3],
45    boundaries: [Boundary; 3],
46    cells: Vec<u32>,
47}
48
49impl VoxelGrid {
50    /// Build a grid of `dims` voxels of edge lengths `spacing_m`, where voxel
51    /// `i + nx (j + ny k)` has material `cells[...]` (an index into
52    /// `materials`), with `boundaries` per axis `[x, y, z]`.
53    ///
54    /// Every dimension must be nonzero, every spacing finite and positive
55    /// with a finite extent `n h`, the voxel count must not overflow, and
56    /// every material index must exist.
57    pub fn new(
58        materials: Vec<Material>,
59        dims: [usize; 3],
60        spacing_m: [f64; 3],
61        cells: Vec<u32>,
62        boundaries: [Boundary; 3],
63    ) -> Result<Self, GeometryError> {
64        if materials.is_empty() {
65            return Err(GeometryError::NoMaterials);
66        }
67        let mut extent_m = [0.0; 3];
68        for axis in 0..3 {
69            if dims[axis] == 0 {
70                return Err(GeometryError::ZeroDimension { axis });
71            }
72            let h = spacing_m[axis];
73            let l = dims[axis] as f64 * h;
74            if !(h.is_finite() && h > 0.0 && l.is_finite()) {
75                return Err(GeometryError::InvalidSpacing { axis, spacing_m: h });
76            }
77            extent_m[axis] = l;
78        }
79        let expected = dims[0]
80            .checked_mul(dims[1])
81            .and_then(|v| v.checked_mul(dims[2]))
82            .ok_or(GeometryError::TooManyVoxels { dims })?;
83        if cells.len() != expected {
84            return Err(GeometryError::CellCountMismatch {
85                expected,
86                actual: cells.len(),
87            });
88        }
89        if let Some((cell, &index)) = cells
90            .iter()
91            .enumerate()
92            .find(|&(_, &m)| m as usize >= materials.len())
93        {
94            return Err(GeometryError::MaterialIndexOutOfRange {
95                cell,
96                index,
97                n_materials: materials.len(),
98            });
99        }
100        Ok(Self {
101            materials,
102            dims,
103            spacing_m,
104            extent_m,
105            boundaries,
106            cells,
107        })
108    }
109
110    /// Voxels per axis `[nx, ny, nz]`.
111    pub fn dims(&self) -> [usize; 3] {
112        self.dims
113    }
114
115    /// Voxel edge lengths `[hx, hy, hz]`, m.
116    pub fn spacing_m(&self) -> [f64; 3] {
117        self.spacing_m
118    }
119
120    /// Box size `[Lx, Ly, Lz]` = `n h`, m.
121    pub fn extent_m(&self) -> [f64; 3] {
122        self.extent_m
123    }
124
125    /// Boundary behaviour of axis `axis` (0 = x, 1 = y, 2 = z).
126    pub fn boundary(&self, axis: usize) -> Boundary {
127        self.boundaries[axis]
128    }
129
130    /// Flat region index of voxel `(i, j, k)`, or `None` out of range.
131    pub fn flat_index(&self, ijk: [usize; 3]) -> Option<usize> {
132        (ijk[0] < self.dims[0] && ijk[1] < self.dims[1] && ijk[2] < self.dims[2])
133            .then(|| ijk[0] + self.dims[0] * (ijk[1] + self.dims[1] * ijk[2]))
134    }
135
136    /// Voxel `(i, j, k)` of a flat region index.
137    pub fn voxel_of(&self, region: usize) -> [usize; 3] {
138        let [nx, ny, _] = self.dims;
139        [region % nx, (region / nx) % ny, region / (nx * ny)]
140    }
141
142    /// Material index of voxel `(i, j, k)`, or `None` out of range.
143    pub fn cell_material(&self, ijk: [usize; 3]) -> Option<usize> {
144        self.flat_index(ijk).map(|r| self.cells[r] as usize)
145    }
146
147    fn flat_unchecked(&self, ijk: [usize; 3]) -> usize {
148        ijk[0] + self.dims[0] * (ijk[1] + self.dims[1] * ijk[2])
149    }
150}
151
152impl Geometry for VoxelGrid {
153    fn n_materials(&self) -> usize {
154        self.materials.len()
155    }
156
157    fn material(&self, index: usize) -> &Material {
158        &self.materials[index]
159    }
160
161    fn n_regions(&self) -> usize {
162        self.cells.len()
163    }
164
165    fn material_index(&self, region: usize) -> usize {
166        self.cells[region] as usize
167    }
168
169    fn locate(&self, pos: [f64; 3]) -> Option<usize> {
170        let mut ijk = [0usize; 3];
171        for a in 0..3 {
172            let x = pos[a];
173            let l = self.extent_m[a];
174            if !(x >= 0.0 && x <= l) {
175                return None;
176            }
177            ijk[a] = if x == l {
178                // The upper face belongs to the periodic image (voxel 0).
179                if self.boundaries[a] == Boundary::Vacuum {
180                    return None;
181                }
182                0
183            } else {
184                // `x / h` can round up to `n` just below the face.
185                ((x / self.spacing_m[a]) as usize).min(self.dims[a] - 1)
186            };
187        }
188        Some(self.flat_unchecked(ijk))
189    }
190
191    fn entry_point(&self) -> [f64; 3] {
192        [0.0, 0.5 * self.extent_m[1], 0.5 * self.extent_m[2]]
193    }
194
195    fn exit(&self, region: usize, pos: [f64; 3], dir: [f64; 3], limit: f64) -> Option<Exit> {
196        match self.flight(region, pos, dir, limit) {
197            Flight::Event(e) => Some(e),
198            Flight::Clear { .. } => None,
199        }
200    }
201
202    /// The traversal behind [`Geometry::exit`]. Without an event, the region
203    /// is the voxel the DDA has reached at the limit: same-material faces up
204    /// to and including the limit have been crossed, so a flight ending
205    /// exactly on a face is in the voxel beyond it along `dir`. A periodic
206    /// wrap is always an event, so a clear flight never wraps.
207    fn flight(&self, region: usize, pos: [f64; 3], dir: [f64; 3], limit: f64) -> Flight {
208        if !dir.iter().all(|d| d.is_finite()) || dir == [0.0; 3] {
209            return Flight::Clear { region };
210        }
211        let home = self.cells[region];
212        let mut cell = self.voxel_of(region);
213        let mut step = [0i64; 3];
214        for a in 0..3 {
215            step[a] = if dir[a] > 0.0 {
216                1
217            } else if dir[a] < 0.0 {
218                -1
219            } else {
220                0
221            };
222        }
223        // Face crossed next on each axis, as a coordinate: the upper face of
224        // the voxel when moving up, the lower face when moving down.
225        let face_coord = |a: usize, c: usize| -> f64 {
226            let k = if step[a] > 0 { c + 1 } else { c };
227            k as f64 * self.spacing_m[a]
228        };
229        let t_of = |a: usize, c: usize| -> f64 {
230            if step[a] == 0 {
231                f64::INFINITY
232            } else {
233                // A start slightly past the face (rounding) crosses at once.
234                ((face_coord(a, c) - pos[a]) / dir[a]).max(0.0)
235            }
236        };
237        let mut t_next = [t_of(0, cell[0]), t_of(1, cell[1]), t_of(2, cell[2])];
238        loop {
239            let t = t_next[0].min(t_next[1]).min(t_next[2]);
240            if t.is_nan() || t > limit {
241                return Flight::Clear {
242                    region: self.flat_unchecked(cell),
243                };
244            }
245            let mut at = [
246                pos[0] + t * dir[0],
247                pos[1] + t * dir[1],
248                pos[2] + t * dir[2],
249            ];
250            let mut wrapped = false;
251            let mut escape: Option<(usize, bool)> = None;
252            let mut new_pos = at;
253            for a in 0..3 {
254                // Tied axes are crossed together (exact equality).
255                if step[a] == 0 || t_next[a] != t {
256                    continue;
257                }
258                at[a] = face_coord(a, cell[a]);
259                new_pos[a] = at[a];
260                let n = self.dims[a];
261                if step[a] > 0 {
262                    if cell[a] + 1 == n {
263                        match self.boundaries[a] {
264                            Boundary::Vacuum => {
265                                escape.get_or_insert((a, true));
266                            }
267                            Boundary::Periodic => {
268                                wrapped = true;
269                                cell[a] = 0;
270                                new_pos[a] = 0.0;
271                            }
272                        }
273                    } else {
274                        cell[a] += 1;
275                    }
276                } else if cell[a] == 0 {
277                    match self.boundaries[a] {
278                        Boundary::Vacuum => {
279                            escape.get_or_insert((a, false));
280                        }
281                        Boundary::Periodic => {
282                            wrapped = true;
283                            cell[a] = n - 1;
284                            new_pos[a] = self.extent_m[a];
285                        }
286                    }
287                } else {
288                    cell[a] -= 1;
289                }
290                t_next[a] = t_of(a, cell[a]);
291            }
292            if let Some((axis, positive)) = escape {
293                let face = match (axis, positive) {
294                    (0, false) => Face::Front,
295                    (0, true) => Face::Back,
296                    _ => Face::Side,
297                };
298                let mut normal = [0.0; 3];
299                normal[axis] = if positive { 1.0 } else { -1.0 };
300                return Flight::Event(Exit {
301                    distance: t,
302                    at,
303                    outcome: ExitOutcome::Escape { face, normal },
304                });
305            }
306            let next = self.flat_unchecked(cell);
307            if wrapped || self.cells[next] != home {
308                // A periodic image is the same voxel set: ownership puts the
309                // wrapped coordinate in the voxel just entered.
310                return Flight::Event(Exit {
311                    distance: t,
312                    at,
313                    outcome: ExitOutcome::Enter {
314                        region: next,
315                        pos: new_pos,
316                    },
317                });
318            }
319        }
320    }
321
322    fn in_vacuum(&self, pos: [f64; 3]) -> bool {
323        (0..3).any(|a| {
324            self.boundaries[a] == Boundary::Vacuum && (pos[a] < 0.0 || pos[a] > self.extent_m[a])
325        })
326    }
327}
328
329#[cfg(test)]
330mod tests {
331    use super::super::test_util::{c, enter, si, unit};
332    use super::super::Stack;
333    use super::*;
334
335    const V: Boundary = Boundary::Vacuum;
336    const P: Boundary = Boundary::Periodic;
337
338    /// `dims` voxels of unit spacing, material `f(i, j, k)` (0 or 1).
339    fn grid(
340        dims: [usize; 3],
341        b: [Boundary; 3],
342        f: impl Fn(usize, usize, usize) -> u32,
343    ) -> VoxelGrid {
344        let mut cells = Vec::new();
345        for k in 0..dims[2] {
346            for j in 0..dims[1] {
347                for i in 0..dims[0] {
348                    cells.push(f(i, j, k));
349                }
350            }
351        }
352        VoxelGrid::new(vec![si(), c()], dims, [1.0; 3], cells, b).unwrap()
353    }
354
355    #[test]
356    fn voxel_validation() {
357        let m = || vec![si()];
358        let ok = |dims, sp: [f64; 3], cells: Vec<u32>| VoxelGrid::new(m(), dims, sp, cells, [V; 3]);
359        assert!(ok([1, 1, 1], [1.0; 3], vec![0]).is_ok());
360        assert_eq!(
361            VoxelGrid::new(vec![], [1; 3], [1.0; 3], vec![0], [V; 3]),
362            Err(GeometryError::NoMaterials)
363        );
364        assert_eq!(
365            ok([1, 0, 1], [1.0; 3], vec![]),
366            Err(GeometryError::ZeroDimension { axis: 1 })
367        );
368        for h in [0.0, -1.0, f64::NAN, f64::INFINITY, f64::MAX] {
369            let r = ok([2, 1, 1], [1.0, 1.0, h], vec![0; 2]);
370            let r = if h == f64::MAX {
371                ok([2, 1, 1], [h, 1.0, 1.0], vec![0; 2])
372            } else {
373                r
374            };
375            assert!(
376                matches!(r, Err(GeometryError::InvalidSpacing { .. })),
377                "{h}: {r:?}"
378            );
379        }
380        assert!(matches!(
381            ok([usize::MAX, 2, 1], [1.0; 3], vec![]),
382            Err(GeometryError::TooManyVoxels { .. })
383        ));
384        assert_eq!(
385            ok([2, 2, 1], [1.0; 3], vec![0; 3]),
386            Err(GeometryError::CellCountMismatch {
387                expected: 4,
388                actual: 3
389            })
390        );
391        assert!(matches!(
392            ok([2, 1, 1], [1.0; 3], vec![0, 1]),
393            Err(GeometryError::MaterialIndexOutOfRange {
394                cell: 1,
395                index: 1,
396                n_materials: 1
397            })
398        ));
399    }
400
401    #[test]
402    fn flattening_and_locate() {
403        let g = grid([3, 2, 2], [V, P, V], |i, j, k| {
404            u32::from(i == 2 && j == 1 && k == 1)
405        });
406        // x fastest.
407        assert_eq!(g.flat_index([2, 1, 1]), Some(11));
408        assert_eq!(g.voxel_of(11), [2, 1, 1]);
409        assert_eq!(g.flat_index([3, 0, 0]), None);
410        assert_eq!(g.cell_material([2, 1, 1]), Some(1));
411        assert_eq!(g.material_index(11), 1);
412        assert_eq!(g.n_regions(), 12);
413        assert_eq!(g.locate([2.5, 1.5, 1.5]), Some(11));
414        // Interior faces belong to the higher voxel.
415        assert_eq!(g.locate([1.0, 0.0, 0.0]), Some(1));
416        // Upper face: vacuum is outside, a periodic axis wraps to voxel 0.
417        assert_eq!(g.locate([3.0, 0.5, 0.5]), None);
418        assert_eq!(g.locate([0.5, 2.0, 0.5]), Some(0));
419        assert_eq!(g.locate([-1e-12, 0.5, 0.5]), None);
420        assert_eq!(g.locate([0.5, 0.5, f64::NAN]), None);
421        // Just below the face never rounds into a nonexistent voxel.
422        assert_eq!(
423            g.locate([f64::from_bits(3.0f64.to_bits() - 1), 0.5, 0.5]),
424            Some(2)
425        );
426        assert_eq!(g.entry_point(), [0.0, 1.0, 1.0]);
427        assert!(g.in_vacuum([-0.1, 0.0, 0.0]));
428        assert!(g.in_vacuum([0.0, 0.0, 2.1]));
429        assert!(!g.in_vacuum([0.0, -5.0, 0.0])); // periodic axis
430    }
431
432    #[test]
433    fn oblique_ray_stops_at_a_material_change_and_skips_equal_voxels() {
434        // Material 1 for x >= 2; y faces are crossed on the way but are inside
435        // one material, so they are not events.
436        let g = grid([4, 3, 1], [V, V, P], |i, _, _| u32::from(i >= 2));
437        let d = unit([2.0, 0.5, 0.0]);
438        let start = [0.5, 0.2, 0.5];
439        let e = g.exit(0, start, d, 100.0).unwrap();
440        let t = 1.5 / d[0];
441        assert!((e.distance - t).abs() < 1e-14 * t);
442        assert_eq!(e.at[0], 2.0);
443        assert!((e.at[1] - (0.2 + t * d[1])).abs() < 1e-14);
444        let (r, pos) = enter(&e);
445        assert_eq!(r, g.flat_index([2, 0, 0]).unwrap());
446        assert_eq!(pos, e.at);
447        // Not reached within a shorter limit; reached at exactly the limit.
448        assert!(g.exit(0, start, d, 0.99 * t).is_none());
449        assert!(g.exit(0, start, d, e.distance).is_some());
450    }
451
452    #[test]
453    fn oblique_ray_changes_material_across_a_lateral_face() {
454        let g = grid([3, 3, 1], [V, V, P], |_, j, _| u32::from(j >= 1));
455        let d = unit([1.0, 2.0, 0.0]);
456        let e = g.exit(0, [0.5, 0.5, 0.5], d, 10.0).unwrap();
457        assert!((e.distance - 0.5 / d[1]).abs() < 1e-14);
458        assert_eq!(e.at[1], 1.0);
459        assert!((e.at[0] - (0.5 + e.distance * d[0])).abs() < 1e-14);
460        assert_eq!(enter(&e).0, g.flat_index([0, 1, 0]).unwrap());
461    }
462
463    #[test]
464    fn negative_direction_and_starting_on_a_face() {
465        let g = grid([4, 1, 1], [V, P, P], |i, _, _| u32::from(i >= 2));
466        // From voxel 3 going back: the face at x = 2.
467        let e = g.exit(3, [3.5, 0.5, 0.5], [-1.0, 0.0, 0.0], 10.0).unwrap();
468        assert_eq!(e.distance, 1.5);
469        assert_eq!(enter(&e), (1, [2.0, 0.5, 0.5]));
470        // Starting on the face x = 2 (owned by voxel 2): moving down crosses
471        // at once, moving up does not see that face.
472        assert_eq!(g.locate([2.0, 0.5, 0.5]), Some(2));
473        let e = g.exit(2, [2.0, 0.5, 0.5], [-1.0, 0.0, 0.0], 10.0).unwrap();
474        assert_eq!(e.distance, 0.0);
475        assert_eq!(enter(&e).0, 1);
476        let e = g.exit(2, [2.0, 0.5, 0.5], [1.0, 0.0, 0.0], 10.0).unwrap();
477        // Same material up to the back face: vacuum escape at x = 4.
478        assert_eq!(e.distance, 2.0);
479        assert_eq!(
480            e.outcome,
481            ExitOutcome::Escape {
482                face: Face::Back,
483                normal: [1.0, 0.0, 0.0]
484            }
485        );
486    }
487
488    #[test]
489    fn edge_and_corner_ties_cross_together() {
490        // Only the diagonal voxel (1,1,0) differs: a 45 degree ray from the
491        // centre of (0,0,0) hits the corner at (1,1) and enters it directly.
492        let g = grid([3, 3, 1], [V, V, P], |i, j, _| u32::from(i == 1 && j == 1));
493        let d = unit([1.0, 1.0, 0.0]);
494        let e = g.exit(0, [0.5, 0.5, 0.5], d, 10.0).unwrap();
495        assert!((e.distance - 0.5 * 2f64.sqrt()).abs() < 1e-15);
496        assert_eq!(e.at[..2], [1.0, 1.0]);
497        assert_eq!(enter(&e).0, g.flat_index([1, 1, 0]).unwrap());
498
499        // Only the two voxels the ray merely touches differ: it passes the
500        // corner without any event and leaves through the far corner, which
501        // is reported as the back face (lowest axis first).
502        let g = grid([3, 3, 1], [V, V, P], |i, j, _| u32::from(i != j));
503        let e = g.exit(0, [0.5, 0.5, 0.5], d, 10.0).unwrap();
504        assert!((e.distance - 2.5 * 2f64.sqrt()).abs() < 1e-14);
505        assert_eq!(
506            e.outcome,
507            ExitOutcome::Escape {
508                face: Face::Back,
509                normal: [1.0, 0.0, 0.0]
510            }
511        );
512        assert_eq!(e.at[..2], [3.0, 3.0]);
513    }
514
515    #[test]
516    fn side_faces_are_reported_as_side() {
517        let g = grid([2, 2, 2], [P, V, V], |_, _, _| 0);
518        let e = g.exit(0, [0.5, 0.5, 0.5], [0.0, -1.0, 0.0], 10.0).unwrap();
519        assert_eq!(
520            e.outcome,
521            ExitOutcome::Escape {
522                face: Face::Side,
523                normal: [0.0, -1.0, 0.0]
524            }
525        );
526        assert_eq!(e.distance, 0.5);
527        let e = g.exit(0, [0.5, 0.5, 0.5], [0.0, 0.0, 1.0], 10.0).unwrap();
528        assert_eq!(e.distance, 1.5);
529        assert!(matches!(
530            e.outcome,
531            ExitOutcome::Escape {
532                face: Face::Side,
533                normal: [0.0, 0.0, 1.0]
534            }
535        ));
536    }
537
538    #[test]
539    fn periodic_wrap_is_an_event_that_keeps_the_path_length() {
540        let g = grid([2, 2, 1], [P, P, P], |_, _, _| 0);
541        let e = g.exit(1, [1.5, 0.5, 0.5], [1.0, 0.0, 0.0], 10.0).unwrap();
542        assert_eq!(e.distance, 0.5);
543        assert_eq!(e.at, [2.0, 0.5, 0.5]);
544        assert_eq!(enter(&e), (0, [0.0, 0.5, 0.5]));
545        let e = g.exit(0, [0.5, 0.5, 0.5], [-1.0, 0.0, 0.0], 10.0).unwrap();
546        assert_eq!(e.distance, 0.5);
547        assert_eq!(enter(&e), (1, [2.0, 0.5, 0.5]));
548        // Oblique, wrapping on y after crossing (same-material) x faces.
549        let d = unit([1.0, 3.0, 0.0]);
550        let e = g.exit(0, [0.1, 0.5, 0.5], d, 10.0).unwrap();
551        assert!((e.distance - 1.5 / d[1]).abs() < 1e-14);
552        let (r, pos) = enter(&e);
553        assert_eq!(pos[1], 0.0);
554        assert!((pos[0] - (0.1 + e.distance * d[0])).abs() < 1e-14);
555        assert_eq!(r, g.flat_index([0, 0, 0]).unwrap());
556        // A single-voxel periodic axis wraps into itself.
557        let g1 = VoxelGrid::new(vec![si()], [1, 1, 1], [2.0; 3], vec![0], [P; 3]).unwrap();
558        let e = g1.exit(0, [0.5, 1.0, 1.0], [-1.0, 0.0, 0.0], 10.0).unwrap();
559        assert_eq!(e.distance, 0.5);
560        assert_eq!(enter(&e), (0, [2.0, 1.0, 1.0]));
561    }
562
563    #[test]
564    fn clear_flight_reports_the_region_it_ends_in() {
565        // Regression: [2, 2, 1] voxels, material 1 only in voxel (1, 1, 0).
566        let g = grid([2, 2, 1], [V; 3], |i, j, _| u32::from(i == 1 && j == 1));
567        let (x, y) = ([1.0, 0.0, 0.0], [0.0, 1.0, 0.0]);
568        // One unit along +x crosses only a same-material face: no event, but
569        // the flight ends in voxel (1, 0, 0).
570        let f = g.flight(0, [0.5; 3], x, 1.0);
571        assert_eq!(f, Flight::Clear { region: 1 });
572        assert!(g.exit(0, [0.5; 3], x, 1.0).is_none());
573        // Turning to +y from there meets the material interface at y = 1.
574        let e = g.exit(1, [1.5, 0.5, 0.5], y, 1.0).unwrap();
575        assert_eq!(e.distance, 0.5);
576        assert_eq!(enter(&e), (3, [1.5, 1.0, 0.5]));
577        // The stale starting region would miss it (x-cell 0 has no y change).
578        assert!(g.exit(0, [1.5, 0.5, 0.5], y, 1.0).is_none());
579
580        // Directional ownership: a flight ending exactly on a face is in the
581        // voxel beyond it along the direction, as for an `Enter`.
582        assert_eq!(g.flight(0, [0.5; 3], x, 0.5), Flight::Clear { region: 1 });
583        let back = [-1.0, 0.0, 0.0];
584        assert_eq!(
585            g.flight(1, [1.5, 0.5, 0.5], back, 0.5),
586            Flight::Clear { region: 0 }
587        );
588        // Short of the face: unchanged. Degenerate direction: unchanged.
589        assert_eq!(g.flight(0, [0.5; 3], x, 0.25), Flight::Clear { region: 0 });
590        assert_eq!(
591            g.flight(1, [1.5; 3], [0.0; 3], 9.0),
592            Flight::Clear { region: 1 }
593        );
594
595        // Periodic: after a wrap the particle sits on the upper face in the
596        // last voxel (`locate` would give voxel 0); a clear flight from there
597        // keeps that ownership, and reaching the wrap is always an event.
598        let p = grid([2, 1, 1], [P, P, P], |_, _, _| 0);
599        assert_eq!(p.locate([2.0, 0.5, 0.5]), Some(0));
600        assert_eq!(
601            p.flight(1, [2.0, 0.5, 0.5], back, 0.25),
602            Flight::Clear { region: 1 }
603        );
604        assert!(matches!(
605            p.flight(1, [1.5, 0.5, 0.5], x, 0.5),
606            Flight::Event(Exit {
607                outcome: ExitOutcome::Enter { region: 0, .. },
608                ..
609            })
610        ));
611
612        // A stack keeps its region: leaving a layer is always an event.
613        let s = Stack::new(vec![(si(), 1.0), (c(), 2.0)], None).unwrap();
614        assert_eq!(s.flight(0, [0.5; 3], x, 0.25), Flight::Clear { region: 0 });
615        assert!(matches!(s.flight(0, [0.5; 3], x, 0.5), Flight::Event(_)));
616    }
617
618    #[test]
619    fn degenerate_directions_have_no_events() {
620        let g = grid([2, 1, 1], [V, V, V], |_, _, _| 0);
621        for d in [[0.0; 3], [f64::NAN, 0.0, 0.0], [f64::INFINITY, 0.0, 0.0]] {
622            assert!(g.exit(0, [0.5; 3], d, f64::INFINITY).is_none());
623        }
624    }
625
626    #[test]
627    fn long_rays_do_not_drift() {
628        // 1000 voxels along x with a 1e-3 tilt: the crossing of the last
629        // x-face is computed from the face coordinate, not accumulated.
630        let n = 1000;
631        let g = VoxelGrid::new(
632            vec![si()],
633            [n, 1, 1],
634            [0.1, 1e3, 1e3],
635            vec![0; n],
636            [V, P, P],
637        )
638        .unwrap();
639        let d = unit([1.0, 1e-3, 0.0]);
640        let e = g.exit(0, [0.05, 0.05, 0.05], d, 1e9).unwrap();
641        assert_eq!(e.at[0], 100.0);
642        assert!((e.distance - (100.0 - 0.05) / d[0]).abs() < 1e-12);
643    }
644}
lindhard/geometry/voxel.rs

lindhard/geometry/
voxel.rs
