Skip to main content

proof_engine/solver/
boundary.rs

1//! Boundary condition handling — Dirichlet, Neumann, periodic, absorbing.
2
3use super::pde::ScalarField2D;
4
5/// Boundary condition type.
6#[derive(Debug, Clone)]
7pub enum BoundaryType {
8    /// Fixed value at boundary.
9    Dirichlet(f64),
10    /// Fixed gradient (normal derivative) at boundary.
11    Neumann(f64),
12    /// Periodic: opposite edges wrap around.
13    Periodic,
14    /// Absorbing: values decay at boundary (open boundary).
15    Absorbing { decay: f64 },
16}
17
18/// Boundary condition for a 2D scalar field.
19#[derive(Debug, Clone)]
20pub struct BoundaryCondition {
21    pub bc_type: BoundaryType,
22}
23
24impl BoundaryCondition {
25    pub fn new(bc_type: BoundaryType) -> Self { Self { bc_type } }
26    pub fn dirichlet(val: f64) -> Self { Self::new(BoundaryType::Dirichlet(val)) }
27    pub fn neumann(grad: f64) -> Self { Self::new(BoundaryType::Neumann(grad)) }
28    pub fn periodic() -> Self { Self::new(BoundaryType::Periodic) }
29    pub fn absorbing(decay: f64) -> Self { Self::new(BoundaryType::Absorbing { decay }) }
30
31    /// Apply boundary conditions to a field.
32    pub fn apply(&self, field: &mut ScalarField2D) {
33        let w = field.width;
34        let h = field.height;
35
36        match &self.bc_type {
37            BoundaryType::Dirichlet(val) => {
38                for x in 0..w { field.set(x, 0, *val); field.set(x, h - 1, *val); }
39                for y in 0..h { field.set(0, y, *val); field.set(w - 1, y, *val); }
40            }
41            BoundaryType::Neumann(grad) => {
42                let dx = field.dx;
43                for x in 1..w - 1 {
44                    field.data[x] = field.data[w + x] - grad * dx;
45                    field.data[(h - 1) * w + x] = field.data[(h - 2) * w + x] + grad * dx;
46                }
47                for y in 1..h - 1 {
48                    field.data[y * w] = field.data[y * w + 1] - grad * dx;
49                    field.data[y * w + w - 1] = field.data[y * w + w - 2] + grad * dx;
50                }
51            }
52            BoundaryType::Periodic => {
53                for x in 0..w {
54                    field.data[x] = field.data[(h - 2) * w + x];
55                    field.data[(h - 1) * w + x] = field.data[w + x];
56                }
57                for y in 0..h {
58                    field.data[y * w] = field.data[y * w + w - 2];
59                    field.data[y * w + w - 1] = field.data[y * w + 1];
60                }
61            }
62            BoundaryType::Absorbing { decay } => {
63                for x in 0..w {
64                    field.data[x] *= 1.0 - decay;
65                    field.data[(h - 1) * w + x] *= 1.0 - decay;
66                }
67                for y in 0..h {
68                    field.data[y * w] *= 1.0 - decay;
69                    field.data[y * w + w - 1] *= 1.0 - decay;
70                }
71            }
72        }
73    }
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79
80    #[test]
81    fn dirichlet_sets_boundary() {
82        let mut field = ScalarField2D::new(5, 5, 1.0, 1.0);
83        field.fill(10.0);
84        BoundaryCondition::dirichlet(0.0).apply(&mut field);
85        assert_eq!(field.get(0, 0), 0.0);
86        assert_eq!(field.get(2, 2), 10.0); // interior unchanged
87    }
88
89    #[test]
90    fn periodic_wraps() {
91        let mut field = ScalarField2D::new(5, 5, 1.0, 1.0);
92        field.set(1, 1, 42.0);
93        field.set(3, 1, 99.0);
94        BoundaryCondition::periodic().apply(&mut field);
95        assert_eq!(field.get(0, 1), field.get(3, 1));
96    }
97}