neco_contact/
lib.rs

1use neco_array2::Array2;
2
3/// Return the index of the cell nearest to `(tx, ty)` on a 2D grid.
4pub fn find_nearest(x: &Array2<f64>, y: &Array2<f64>, tx: f64, ty: f64) -> (usize, usize) {
5    let (nx, ny) = x.dim();
6    let mut best = (0, 0);
7    let mut best_distance = f64::INFINITY;
8    for i in 0..nx {
9        for j in 0..ny {
10            let distance = (x[[i, j]] - tx).powi(2) + (y[[i, j]] - ty).powi(2);
11            if distance < best_distance {
12                best_distance = distance;
13                best = (i, j);
14            }
15        }
16    }
17    best
18}
19
20/// Build a normalized cosine-taper spatial mask.
21pub fn build_spatial_mask(
22    x: &Array2<f64>,
23    y: &Array2<f64>,
24    hx: f64,
25    hy: f64,
26    width: f64,
27    interior: Option<&Array2<bool>>,
28) -> Array2<f64> {
29    let (nx, ny) = x.dim();
30    let mut mask = Array2::zeros((nx, ny));
31    let mut sum = 0.0;
32    for i in 0..nx {
33        for j in 0..ny {
34            if let Some(active) = interior {
35                if !active[[i, j]] {
36                    continue;
37                }
38            }
39            let distance = ((x[[i, j]] - hx).powi(2) + (y[[i, j]] - hy).powi(2)).sqrt();
40            if distance < width {
41                let value = 0.5 * (1.0 + (std::f64::consts::PI * distance / width).cos());
42                mask[[i, j]] = value;
43                sum += value;
44            }
45        }
46    }
47    if sum > 0.0 {
48        for value in mask.iter_mut() {
49            *value /= sum;
50        }
51    }
52    mask
53}
54
55/// Collect active interior cells that are at least `margin` cells from the border.
56pub fn collect_interior(interior: &Array2<bool>, margin: usize) -> Vec<(usize, usize)> {
57    let (nx, ny) = interior.dim();
58    let mut points = Vec::new();
59    for i in margin..nx.saturating_sub(margin) {
60        for j in margin..ny.saturating_sub(margin) {
61            if interior[[i, j]] {
62                points.push((i, j));
63            }
64        }
65    }
66    points
67}
68
69/// Common Hertz contact update for beater-based excitation.
70#[derive(Debug, Clone)]
71pub struct HertzContact {
72    pub beater_x: f64,
73    pub beater_v: f64,
74    beater_mass: f64,
75    k_hertz: f64,
76    alpha_hertz: f64,
77    contact_ended: bool,
78}
79
80impl HertzContact {
81    pub fn new(beater_mass: f64, k_hertz: f64, alpha_hertz: f64, v0: f64) -> Self {
82        Self {
83            beater_x: 0.0,
84            beater_v: v0,
85            beater_mass,
86            k_hertz,
87            alpha_hertz,
88            contact_ended: false,
89        }
90    }
91
92    /// Advance one step and return the contact force.
93    pub fn step(&mut self, w_surface: f64, dt: f64) -> f64 {
94        if self.contact_ended {
95            return 0.0;
96        }
97        let delta = self.beater_x - w_surface;
98        if delta > 0.0 {
99            let force = self.k_hertz * delta.powf(self.alpha_hertz);
100            let acceleration = -force / self.beater_mass;
101            self.beater_v += acceleration * dt;
102            self.beater_x += self.beater_v * dt;
103            force
104        } else {
105            self.beater_x += self.beater_v * dt;
106            if self.beater_x > 1e-12 && delta < -1e-10 {
107                self.contact_ended = true;
108            }
109            0.0
110        }
111    }
112
113    pub fn energy(&self) -> f64 {
114        0.5 * self.beater_mass * self.beater_v * self.beater_v
115    }
116
117    pub fn contact_ended(&self) -> bool {
118        self.contact_ended
119    }
120
121    pub fn set_contact_ended(&mut self, ended: bool) {
122        self.contact_ended = ended;
123    }
124}
neco_contact/lib.rs

neco_contact/
lib.rs
