axiolid_curve/
pair_section.rs

1//! Where two parametric surfaces meet, carried by nodes on both (ADR 0077).
2//!
3//! Two B-spline surfaces have no implicit equation to read one in the
4//! other's parameters, so their section is carried directly: a chain of
5//! [`PairNode`]s, each a point on both surfaces with its parameters on each,
6//! found to the last bits. Between two nodes the curve is defined, not
7//! interpolated: at local parameter `s` it is the point where both surfaces
8//! meet on the plane across the chord at `P0 + s (P1 - P0)`, the unique
9//! solution of four equations in the four parameters near the chord --
10//! `S1(a) = S2(b)` and the plane -- found by Newton from the nodes'
11//! parameters. Derivatives follow from the same system.
12
13use axiolid_core::{Point2, Point3, Scalar, Vec3};
14
15use crate::implicit::{Carrier, SurfaceJet};
16
17/// A point on both surfaces, with its parameters on each.
18#[derive(Debug, Clone, Copy, PartialEq)]
19#[cfg_attr(
20    feature = "serde",
21    derive(serde::Serialize, serde::Deserialize),
22    serde(deny_unknown_fields)
23)]
24pub struct PairNode {
25    /// The point.
26    pub point: Point3,
27    /// Its parameters on the first surface.
28    pub first: Point2,
29    /// Its parameters on the second surface.
30    pub second: Point2,
31}
32
33/// A stretch of the curve where two surfaces meet. The parameter runs over
34/// `[0, nodes.len() - 1]`, one unit per chord.
35#[derive(Debug, Clone, PartialEq)]
36#[cfg_attr(
37    feature = "serde",
38    derive(serde::Serialize, serde::Deserialize),
39    serde(deny_unknown_fields)
40)]
41pub struct PairSection3 {
42    /// The first surface.
43    pub first: Carrier,
44    /// The second surface.
45    pub second: Carrier,
46    /// Nodes on both, in order along the curve.
47    pub nodes: Vec<PairNode>,
48}
49
50/// Solve the 4x4 system `m x = r` by Gaussian elimination with partial
51/// pivoting; `None` when it is singular.
52#[must_use]
53#[allow(clippy::needless_range_loop)]
54pub fn solve4(mut m: [[Scalar; 4]; 4], mut r: [Scalar; 4]) -> Option<[Scalar; 4]> {
55    for col in 0..4 {
56        let pivot = (col..4).max_by(|&a, &b| m[a][col].abs().total_cmp(&m[b][col].abs()))?;
57        if m[pivot][col] == 0.0 || !m[pivot][col].is_finite() {
58            return None;
59        }
60        m.swap(col, pivot);
61        r.swap(col, pivot);
62        for row in col + 1..4 {
63            let f = m[row][col] / m[col][col];
64            for k in col..4 {
65                m[row][k] -= f * m[col][k];
66            }
67            r[row] -= f * r[col];
68        }
69    }
70    let mut x = [0.0; 4];
71    for row in (0..4).rev() {
72        let mut acc = r[row];
73        for k in row + 1..4 {
74            acc -= m[row][k] * x[k];
75        }
76        x[row] = acc / m[row][row];
77    }
78    x.iter().all(|v| v.is_finite()).then_some(x)
79}
80
81/// The system's matrix at `(a, b)` for chord direction `d`.
82fn matrix(j1: &SurfaceJet, j2: &SurfaceJet, d: Vec3) -> [[Scalar; 4]; 4] {
83    [
84        [j1.u.x, j1.v.x, -j2.u.x, -j2.v.x],
85        [j1.u.y, j1.v.y, -j2.u.y, -j2.v.y],
86        [j1.u.z, j1.v.z, -j2.u.z, -j2.v.z],
87        [d.dot(j1.u), d.dot(j1.v), 0.0, 0.0],
88    ]
89}
90
91impl PairSection3 {
92    /// The parameter range's end, `nodes.len() - 1`.
93    #[must_use]
94    pub fn end(&self) -> Scalar {
95        self.nodes.len().saturating_sub(1) as Scalar
96    }
97
98    fn locate(&self, t: Scalar) -> Option<(usize, Scalar)> {
99        let n = self.nodes.len();
100        if n < 2 || !t.is_finite() {
101            return None;
102        }
103        let slack = 1e-12 * (1.0 + self.end());
104        if t < -slack || t > self.end() + slack {
105            return None;
106        }
107        let i = (t.floor().max(0.0) as usize).min(n - 2);
108        Some((i, (t - i as Scalar).clamp(0.0, 1.0)))
109    }
110
111    /// The parameters on both surfaces at `t`, and the point.
112    #[must_use]
113    pub fn solve(&self, t: Scalar) -> Option<(Point2, Point2, Point3)> {
114        let (i, s) = self.locate(t)?;
115        let (n0, n1) = (self.nodes[i], self.nodes[i + 1]);
116        if s == 0.0 {
117            return Some((n0.first, n0.second, n0.point));
118        }
119        if s == 1.0 {
120            return Some((n1.first, n1.second, n1.point));
121        }
122        let d = n1.point - n0.point;
123        let target = n0.point + d * s;
124        let mut a = n0.first + (n1.first - n0.first) * s;
125        let mut b = n0.second + (n1.second - n0.second) * s;
126        let scale = 1.0 + n0.point.length() + d.length();
127        for _ in 0..50 {
128            let j1 = self.first.jet(a.x, a.y);
129            let j2 = self.second.jet(b.x, b.y);
130            let gap = j1.point - j2.point;
131            let plane = d.dot(j1.point - target);
132            let x = solve4(matrix(&j1, &j2, d), [-gap.x, -gap.y, -gap.z, -plane])?;
133            a += Point2::new(x[0], x[1]);
134            b += Point2::new(x[2], x[3]);
135            let step = x.iter().map(|v| v.abs()).fold(0.0, Scalar::max);
136            if step <= 4.0 * Scalar::EPSILON * (1.0 + a.length() + b.length()) {
137                break;
138            }
139        }
140        let (j1, j2) = (self.first.jet(a.x, a.y), self.second.jet(b.x, b.y));
141        ((j1.point - j2.point).length() <= 1e-9 * scale).then_some((a, b, j1.point))
142    }
143
144    /// The point at `t`.
145    #[must_use]
146    pub fn point(&self, t: Scalar) -> Option<Point3> {
147        self.solve(t).map(|(_, _, p)| p)
148    }
149
150    /// The rates of both surfaces' parameters and of the point at `t`:
151    /// from `S1_a a' - S2_b b' = 0` and `d . S1_a a' = |d|^2`.
152    #[must_use]
153    pub fn rates(&self, t: Scalar) -> Option<(Point2, Point2, Vec3)> {
154        let (i, _) = self.locate(t)?;
155        let d = self.nodes[i + 1].point - self.nodes[i].point;
156        let (a, b, _) = self.solve(t)?;
157        let (j1, j2) = (self.first.jet(a.x, a.y), self.second.jet(b.x, b.y));
158        let x = solve4(matrix(&j1, &j2, d), [0.0, 0.0, 0.0, d.dot(d)])?;
159        let (da, db) = (Point2::new(x[0], x[1]), Point2::new(x[2], x[3]));
160        Some((da, db, j1.u * da.x + j1.v * da.y))
161    }
162
163    /// `dP/dt`.
164    #[must_use]
165    pub fn tangent(&self, t: Scalar) -> Option<Vec3> {
166        self.rates(t).map(|(_, _, d)| d)
167    }
168
169    /// The second rates of both surfaces' parameters and of the point at
170    /// `t`, by central differences of [`Self::rates`] inside the chord.
171    #[must_use]
172    pub fn second_rates(&self, t: Scalar) -> Option<(Point2, Point2, Vec3)> {
173        let (i, s) = self.locate(t)?;
174        let h = 1e-5;
175        let (lo, hi) = ((s - h).max(0.0), (s + h).min(1.0));
176        let base = i as Scalar;
177        let (p, q) = (self.rates(base + lo)?, self.rates(base + hi)?);
178        let w = hi - lo;
179        Some(((q.0 - p.0) / w, (q.1 - p.1) / w, (q.2 - p.2) / w))
180    }
181
182    /// `d2P/dt2`.
183    #[must_use]
184    pub fn bend(&self, t: Scalar) -> Option<Vec3> {
185        self.second_rates(t).map(|(_, _, d)| d)
186    }
187
188    /// Which of the two surfaces `carrier` is: `Some(true)` for the first,
189    /// `Some(false)` for the second.
190    #[must_use]
191    pub fn side(&self, carrier: &Carrier) -> Option<bool> {
192        if *carrier == self.first {
193            Some(true)
194        } else if *carrier == self.second {
195            Some(false)
196        } else {
197            None
198        }
199    }
200
201    /// The stretch from `t0` to `t1` (`t0 < t1`) as a curve of its own over
202    /// `[0, nodes]`: the nodes between, and new end nodes solved at `t0` and
203    /// `t1`.
204    #[must_use]
205    pub fn sub(&self, t0: Scalar, t1: Scalar) -> Option<Self> {
206        let (a0, b0, p0) = self.solve(t0)?;
207        let (a1, b1, p1) = self.solve(t1)?;
208        let mut nodes = vec![PairNode {
209            point: p0,
210            first: a0,
211            second: b0,
212        }];
213        let first = t0.floor() as usize + 1;
214        let last = t1.ceil() as usize;
215        for i in first..last.min(self.nodes.len()) {
216            let t = i as Scalar;
217            if t > t0 + 1e-9 && t < t1 - 1e-9 {
218                nodes.push(self.nodes[i]);
219            }
220        }
221        nodes.push(PairNode {
222            point: p1,
223            first: a1,
224            second: b1,
225        });
226        Some(Self {
227            first: self.first.clone(),
228            second: self.second.clone(),
229            nodes,
230        })
231    }
232
233    /// The same curve run backwards.
234    #[must_use]
235    pub fn reversed(&self) -> Self {
236        let mut nodes = self.nodes.clone();
237        nodes.reverse();
238        Self {
239            first: self.first.clone(),
240            second: self.second.clone(),
241            nodes,
242        }
243    }
244
245    /// The parameter of a point on the curve: the chord nearest it, the
246    /// plane across that chord through the point, then Newton along the
247    /// curve's tangent within the chord.
248    #[must_use]
249    pub fn parameter_of(&self, p: Point3) -> Option<Scalar> {
250        let mut best: Option<(Scalar, usize, Scalar)> = None;
251        for (i, w) in self.nodes.windows(2).enumerate() {
252            let d = w[1].point - w[0].point;
253            let l2 = d.length_squared();
254            if l2 == 0.0 {
255                continue;
256            }
257            let s = ((p - w[0].point).dot(d) / l2).clamp(0.0, 1.0);
258            let miss = (w[0].point + d * s - p).length();
259            if best.is_none_or(|(m, _, _)| miss < m) {
260                best = Some((miss, i, s));
261            }
262        }
263        let (_, i, mut s) = best?;
264        let base = i as Scalar;
265        for _ in 0..30 {
266            let t = base + s;
267            let (q, d) = (self.point(t)?, self.tangent(t)?);
268            let l2 = d.length_squared();
269            if l2 == 0.0 {
270                break;
271            }
272            let next = (s - (q - p).dot(d) / l2).clamp(0.0, 1.0);
273            let step = (next - s).abs();
274            s = next;
275            if step <= 4.0 * Scalar::EPSILON * (1.0 + base) {
276                break;
277            }
278        }
279        Some(base + s)
280    }
281
282    /// Whether every number is finite.
283    #[must_use]
284    pub fn is_finite(&self) -> bool {
285        self.first.is_finite()
286            && self.second.is_finite()
287            && self
288                .nodes
289                .iter()
290                .all(|n| n.point.is_finite() && n.first.is_finite() && n.second.is_finite())
291    }
292}
axiolid_curve/pair_section.rs

axiolid_curve/
pair_section.rs
