geop_core_geometry/nurb_surface/
join.rs

1//! Two surfaces, one running on from where the other ends, as one surface.
2
3use geop_core_math::{
4    geop_error::{GeopError, GeopResult},
5    vector::Vector,
6};
7
8use super::NurbSurface3D;
9use crate::knot_insertion::is_clamped_at;
10
11impl NurbSurface3D {
12    /// This surface and `next`, which starts along `u` — along `v`, if not
13    /// `along_u` — where this one ends, as one surface: `next`'s knots
14    /// carried on from where this one's end, its homogeneous control points
15    /// scaled so that its first row's weights are this one's last row's —
16    /// the same rational surface — and the two joined by a knot of full
17    /// multiplicity, so the joined surface is exactly the two pieces,
18    /// however they meet. The row where they meet is what both had there:
19    /// two enclosures of one row of points, so what they share. What the
20    /// pieces of a wall swept along a profile are, one face per piece of
21    /// the profile, joined into the wall.
22    ///
23    /// Fails, saying why, for surfaces of other degrees or other knots
24    /// across the join, knot vectors not clamped where they meet, and rows
25    /// that do not meet.
26    pub fn joined(&self, next: &Self, along_u: bool) -> GeopResult<Self> {
27        if !along_u {
28            return Ok(self
29                .transposed()
30                .joined(&next.transposed(), true)?
31                .transposed());
32        }
33        let (a, b) = (self, next);
34        let fail = |why: String| Err(GeopError::new(format!("NurbSurface::joined: {why}")));
35        if a.degree_u != b.degree_u || a.degree_v != b.degree_v {
36            return fail(format!(
37                "degrees {}x{} and {}x{} differ",
38                a.degree_u, a.degree_v, b.degree_u, b.degree_v
39            ));
40        }
41        let across = a.knot_vector_v.len() == b.knot_vector_v.len()
42            && (a.knot_vector_v.iter().zip(&b.knot_vector_v)).all(|(x, y)| x.could_be_equal(*y));
43        if !across {
44            return fail(format!(
45                "the knots across the join differ: {:?} and {:?}",
46                a.knot_vector_v, b.knot_vector_v
47            ));
48        }
49        let p = a.degree_u;
50        if !is_clamped_at(&a.knot_vector_u, p, false) || !is_clamped_at(&b.knot_vector_u, p, true) {
51            return fail(format!(
52                "the knot vectors {:?} and {:?} are not clamped where the surfaces meet",
53                a.knot_vector_u, b.knot_vector_u
54            ));
55        }
56        let nv = a.num_v;
57        let last = &a.control_points[(a.num_u - 1) * nv..];
58        let first = &b.control_points[..nv];
59        let scale = last[0][3].div(first[0][3])?;
60        let scaled: Vec<Vector<4>> = (b.control_points.iter())
61            .map(|cp| cp.prod_scalar(scale))
62            .collect();
63        if !(last.iter().zip(&scaled[..nv])).all(|(p, q)| p.could_be_equal(q)) {
64            return fail(format!(
65                "the surfaces do not meet: one ends at {last:?}, the other starts at {:?}",
66                &scaled[..nv]
67            ));
68        }
69        let meet: Vec<Vector<4>> = (last.iter().zip(&scaled[..nv]))
70            .map(|(p, q)| Vector::from_array(std::array::from_fn(|k| p[k].intersect(q[k]))))
71            .collect();
72        let (_, a1) = a.domain_u();
73        let (b0, _) = b.domain_u();
74        let shift = a1.sub(b0);
75        let mut knots: Vec<_> = a.knot_vector_u[..a.knot_vector_u.len() - 1].to_vec();
76        knots.extend(b.knot_vector_u[p + 1..].iter().map(|k| k.add(shift)));
77        let mut points = a.control_points[..(a.num_u - 1) * nv].to_vec();
78        points.extend(meet);
79        points.extend_from_slice(&scaled[nv..]);
80        let across: Vec<_> = (a.knot_vector_v.iter().zip(&b.knot_vector_v))
81            .map(|(x, y)| x.union(*y))
82            .collect();
83        NurbSurface3D::try_new(p, a.degree_v, points, knots, across)
84    }
85}
86
87#[cfg(test)]
88mod tests {
89    use geop_core_math::{scalars::Scalar, vector::Vector3};
90
91    use crate::nurb_surface::NurbSurface3D;
92
93    fn p(x: f64, y: f64, z: f64) -> Vector3 {
94        Vector3::from_array([x, y, z].map(Scalar::from_f64))
95    }
96
97    /// Two squares side by side, joined along `u` and along `v`: one
98    /// surface through both, each half where its piece was.
99    #[test]
100    fn joins_two_squares_along_either_direction() {
101        let square = |x: f64| {
102            NurbSurface3D::bilinear(
103                p(x, 0., 0.),
104                p(x + 1., 0., 0.),
105                p(x + 1., 1., 0.),
106                p(x, 1., 0.),
107            )
108            .unwrap()
109        };
110        let joined = square(0.0).joined(&square(1.0), true).unwrap();
111        let (u0, u1) = joined.domain_u();
112        assert_eq!(u1.sub(u0).to_f64(), 2.0);
113        let at = joined
114            .evaluate(Scalar::from_f64(1.5), Scalar::from_f64(0.5))
115            .unwrap();
116        assert!((at[0].to_f64() - 1.5).abs() < 1e-12 && (at[1].to_f64() - 0.5).abs() < 1e-12);
117        let turned = |s: &NurbSurface3D| s.transposed();
118        let joined = turned(&square(0.0))
119            .joined(&turned(&square(1.0)), false)
120            .unwrap();
121        let at = joined
122            .evaluate(Scalar::from_f64(0.5), Scalar::from_f64(1.5))
123            .unwrap();
124        assert!((at[0].to_f64() - 1.5).abs() < 1e-12, "{at:?}");
125        assert!(square(0.0).joined(&square(2.0), true).is_err(), "apart");
126    }
127}
geop_core_geometry/nurb_surface/join.rs

geop_core_geometry/nurb_surface/
join.rs
