brep_kernel/props/mass_properties/
integration.rs

1use super::*;
2
3pub(super) fn interior_knots(knots: &[f64], degree: usize) -> Vec<f64> {
4    let start = knots[degree];
5    let end = knots[knots.len() - 1 - degree];
6    let mut result = Vec::new();
7    for &knot in knots {
8        if knot <= start + 1e-12 || knot >= end - 1e-12 {
9            continue;
10        }
11        if result
12            .last()
13            .is_none_or(|previous: &f64| (knot - *previous).abs() > 1e-12)
14        {
15            result.push(knot);
16        }
17    }
18    result
19}
20
21pub(super) fn curve_breaks(curve: &NurbsCurve) -> Result<Vec<f64>, String> {
22    let [start, end] = curve.domain()?;
23    let mut result = vec![start];
24    result.extend(interior_knots(&curve.knots, curve.degree));
25    result.push(end);
26    Ok(result)
27}
28
29pub fn parameter_space_area(face: &FaceRecord) -> Result<f64, String> {
30    let mut area = 0.0;
31    for loop_record in &face.loops {
32        for coedge in &loop_record.coedges {
33            for pair in curve_breaks(&coedge.pcurve)?.windows(2) {
34                let half = (pair[1] - pair[0]) * 0.5;
35                let middle = (pair[1] + pair[0]) * 0.5;
36                for index in 0..GAUSS_X.len() {
37                    let parameter = middle + half * GAUSS_X[index];
38                    let (point, tangent) = coedge.pcurve.deriv1(parameter)?;
39                    area +=
40                        GAUSS_W[index] * half * 0.5 * (point.x * tangent.y - point.y * tangent.x);
41                }
42            }
43        }
44    }
45    Ok(area)
46}
47
48pub(super) fn is_affine(surface: &NurbsSurface) -> Result<bool, String> {
49    if surface.degree_u != 1
50        || surface.degree_v != 1
51        || surface.control_points.len() != 2
52        || surface.control_points[0].len() != 2
53    {
54        return Ok(false);
55    }
56    let points = &surface.control_points;
57    let weight = points[0][0].w;
58    if [points[0][1], points[1][0], points[1][1]]
59        .iter()
60        .any(|point| (point.w - weight).abs() > 1e-12)
61    {
62        return Ok(false);
63    }
64    let p00 = points[0][0].point()?;
65    let p01 = points[0][1].point()?;
66    let p10 = points[1][0].point()?;
67    let p11 = points[1][1].point()?;
68    let gap = p11.sub(p10).sub(p01.sub(p00)).length();
69    Ok(gap <= 1e-9 * (1.0 + p00.length().max(p11.length())))
70}
71
72pub(super) fn surface_breaks(surface: &NurbsSurface) -> Result<(Vec<f64>, Vec<f64>), String> {
73    let ku = crate::KnotVector::new(surface.knots_u.clone(), surface.degree_u)?;
74    let kv = crate::KnotVector::new(surface.knots_v.clone(), surface.degree_v)?;
75    let [u0, u1] = ku.domain();
76    let [v0, v1] = kv.domain();
77    let mut u = vec![u0];
78    u.extend(interior_knots(&surface.knots_u, surface.degree_u));
79    u.push(u1);
80    let mut v = vec![v0];
81    v.extend(interior_knots(&surface.knots_v, surface.degree_v));
82    v.push(v1);
83    Ok((u, v))
84}
85
86pub(super) fn integrand_value(kind: Integrand, point: Vec3, weighted_normal: Vec3) -> f64 {
87    let (x, y, z) = (point.x, point.y, point.z);
88    match kind {
89        Integrand::Area => weighted_normal.length(),
90        Integrand::Volume => point.dot(weighted_normal),
91        Integrand::VolumeAbout(reference) => point.sub(reference).dot(weighted_normal),
92        Integrand::MomentX => 0.5 * x * x * weighted_normal.x,
93        Integrand::MomentY => 0.5 * y * y * weighted_normal.y,
94        Integrand::MomentZ => 0.5 * z * z * weighted_normal.z,
95        Integrand::SecondXX => x * x * x / 3.0 * weighted_normal.x,
96        Integrand::SecondYY => y * y * y / 3.0 * weighted_normal.y,
97        Integrand::SecondZZ => z * z * z / 3.0 * weighted_normal.z,
98        Integrand::ProductXY => 0.5 * x * x * y * weighted_normal.x,
99        Integrand::ProductXZ => 0.5 * x * x * z * weighted_normal.x,
100        Integrand::ProductYZ => 0.5 * y * y * z * weighted_normal.y,
101    }
102}
103
104pub(super) fn evaluate_integrand(face: &FaceRecord, u: f64, v: f64, kind: Integrand) -> Result<f64, String> {
105    let (point, su, sv) = face.surface.deriv1(u, v)?;
106    let sign = if face.same_sense { 1.0 } else { -1.0 };
107    let weighted_normal = su.cross(sv).scale(sign);
108    Ok(integrand_value(kind, point, weighted_normal))
109}
110
111pub(super) fn integrate_untrimmed(face: &FaceRecord, kind: Integrand) -> Result<f64, String> {
112    let (u_breaks, v_breaks) = surface_breaks(&face.surface)?;
113    let mut total = 0.0;
114    for upair in u_breaks.windows(2) {
115        let half_u = (upair[1] - upair[0]) * 0.5;
116        let middle_u = (upair[1] + upair[0]) * 0.5;
117        for vpair in v_breaks.windows(2) {
118            let half_v = (vpair[1] - vpair[0]) * 0.5;
119            let middle_v = (vpair[1] + vpair[0]) * 0.5;
120            for i in 0..GAUSS_X.len() {
121                for j in 0..GAUSS_X.len() {
122                    total += GAUSS_W[i]
123                        * GAUSS_W[j]
124                        * half_u
125                        * half_v
126                        * evaluate_integrand(
127                            face,
128                            middle_u + half_u * GAUSS_X[i],
129                            middle_v + half_v * GAUSS_X[j],
130                            kind,
131                        )?;
132                }
133            }
134        }
135    }
136    Ok(total)
137}
138
139/// A bi-periodic band face (surface closed in both u and v, exactly two loops
140/// each a full-wrap constant-cross-level rim) whose trim is NOT expressed with a
141/// seam ruling — the fillet-torus bands OCC/STEP emit as two rim circles. Such a
142/// face integrates to ZERO on the trimmed path (both loops are degenerate iso
143/// lines in parameter space), so it is handled analytically over the seam-cut
144/// rectangle instead. See [`biperiodic_band_range`].
145#[derive(Clone, Copy)]
146pub(super) struct BiBand {
147    /// Which parameter is the periodic (full-wrap) one.
148    p_is_u: bool,
149    /// The two rim cross-levels, ascending (q_lo <= q_hi) in the cross param.
150    q_lo: f64,
151    q_hi: f64,
152    /// True when the material band is the CROSS-SEAM COMPLEMENT of [q_lo, q_hi]
153    /// rather than the between-rims strip (mirrors the tessellator's
154    /// `close_periodic_trim_dir` orientation test).
155    complement: bool,
156}
157
158/// Detect the bi-periodic band configuration of `face` and choose which of the
159/// two regions its rims bound is material. The choice replicates the watertight
160/// tessellator (`watertight_tessellation::close_periodic_trim_dir`): the
161/// material is to the LEFT of the boundary traversal, so the between-rims strip
162/// is kept unless the cross direction wraps AND the rim senses conclusively name
163/// the complement. Returns None for anything that is not a clean two-rim band.
164pub(super) fn biperiodic_band_range(face: &FaceRecord) -> Result<Option<BiBand>, String> {
165    let surface = &face.surface;
166    let (closed_u, closed_v) = surface.closed_directions()?;
167    if !(closed_u && closed_v) {
168        return Ok(None);
169    }
170    let [u0, u1] = surface.domain_u()?;
171    let [v0, v1] = surface.domain_v()?;
172    let u_span = (u1 - u0).abs().max(1e-30);
173    let v_span = (v1 - v0).abs().max(1e-30);
174    // Sample each loop's pcurves (in coedge/traversal order) into (u,v) points.
175    // A DEGENERATE loop — a pole/apex vertex whose whole trace collapses to one
176    // parameter point (a fat fillet torus touches its outer equator at a single
177    // seam point: ABC 00000039 faces 68/70) — is not a band boundary; drop it
178    // and keep the two real rims. Anything else leaves the band ambiguous.
179    let mut loop_points: Vec<Vec<[f64; 2]>> = Vec::with_capacity(2);
180    for loop_record in &face.loops {
181        let mut points = Vec::new();
182        for coedge in &loop_record.coedges {
183            let [d0, d1] = coedge.pcurve.domain()?;
184            let samples = 12;
185            for k in 0..=samples {
186                let t = d0 + (d1 - d0) * k as f64 / samples as f64;
187                let p = coedge.pcurve.evaluate(t)?;
188                points.push([p.x, p.y]);
189            }
190        }
191        if points.len() < 2 {
192            return Ok(None);
193        }
194        let (mut umin, mut umax, mut vmin, mut vmax) = (
195            f64::INFINITY,
196            f64::NEG_INFINITY,
197            f64::INFINITY,
198            f64::NEG_INFINITY,
199        );
200        for pt in &points {
201            umin = umin.min(pt[0]);
202            umax = umax.max(pt[0]);
203            vmin = vmin.min(pt[1]);
204            vmax = vmax.max(pt[1]);
205        }
206        if (umax - umin) <= 1e-3 * u_span && (vmax - vmin) <= 1e-3 * v_span {
207            continue; // degenerate pole/apex loop
208        }
209        loop_points.push(points);
210    }
211    if loop_points.len() != 2 {
212        return Ok(None);
213    }
214    for p_is_u in [true, false] {
215        let (period, _p0, _p1) = if p_is_u {
216            (u1 - u0, u0, u1)
217        } else {
218            (v1 - v0, v0, v1)
219        };
220        let (q_dom_lo, q_dom_hi) = if p_is_u { (v0, v1) } else { (u0, u1) };
221        let q_extent = (q_dom_hi - q_dom_lo).abs().max(1e-30);
222        if !(period > 0.0) {
223            continue;
224        }
225        let coord = |pt: &[f64; 2]| -> (f64, f64) {
226            if p_is_u {
227                (pt[0], pt[1])
228            } else {
229                (pt[1], pt[0])
230            }
231        };
232        // (level, winding-direction) for each rim; None if any loop is not a
233        // clean full-wrap constant-cross-level rim in this p direction.
234        let mut rings: Vec<(f64, i32)> = Vec::with_capacity(2);
235        let mut clean = true;
236        for points in &loop_points {
237            let (mut pmin, mut pmax, mut qmin, mut qmax) = (
238                f64::INFINITY,
239                f64::NEG_INFINITY,
240                f64::INFINITY,
241                f64::NEG_INFINITY,
242            );
243            for pt in points {
244                let (p, q) = coord(pt);
245                pmin = pmin.min(p);
246                pmax = pmax.max(p);
247                qmin = qmin.min(q);
248                qmax = qmax.max(q);
249            }
250            if (pmax - pmin) < 0.6 * period || (qmax - qmin) > 0.05 * q_extent {
251                clean = false;
252                break;
253            }
254            let mut net = 0.0;
255            for pair in points.windows(2) {
256                let mut delta = coord(&pair[1]).0 - coord(&pair[0]).0;
257                if delta > 0.5 * period {
258                    delta -= period;
259                } else if delta < -0.5 * period {
260                    delta += period;
261                }
262                net += delta;
263            }
264            let direction = if net > 0.25 * period {
265                1
266            } else if net < -0.25 * period {
267                -1
268            } else {
269                0
270            };
271            rings.push((0.5 * (qmin + qmax), direction));
272        }
273        if !clean || rings.len() != 2 {
274            continue;
275        }
276        rings.sort_by(|a, b| a.0.total_cmp(&b.0));
277        let (q_lo, lower_dir) = rings[0];
278        let (q_hi, upper_dir) = rings[1];
279        let inconclusive = lower_dir == 0 || upper_dir == 0 || lower_dir == upper_dir;
280        let between_rims_is_ccw_uv = if p_is_u { lower_dir > 0 } else { lower_dir < 0 };
281        let complement = !inconclusive && between_rims_is_ccw_uv != face.same_sense;
282        return Ok(Some(BiBand {
283            p_is_u,
284            q_lo,
285            q_hi,
286            complement,
287        }));
288    }
289    Ok(None)
290}
291
292/// Gauss–Legendre integral of `kind` over the axis-aligned parameter rectangle
293/// [u_lo,u_hi]×[v_lo,v_hi], subdividing at knot spans clipped to the rectangle.
294pub(super) fn integrate_rectangle(
295    face: &FaceRecord,
296    u_lo: f64,
297    u_hi: f64,
298    v_lo: f64,
299    v_hi: f64,
300    kind: Integrand,
301) -> Result<f64, String> {
302    let (u_full, v_full) = surface_breaks(&face.surface)?;
303    let clamp = |breaks: &[f64], lo: f64, hi: f64| -> Vec<f64> {
304        let eps = 1e-9 * (hi - lo).abs().max(1e-30);
305        let mut out = vec![lo];
306        for &b in breaks {
307            if b > lo + eps && b < hi - eps {
308                out.push(b);
309            }
310        }
311        out.push(hi);
312        out
313    };
314    let u_breaks = clamp(&u_full, u_lo, u_hi);
315    let v_breaks = clamp(&v_full, v_lo, v_hi);
316    let mut total = 0.0;
317    for upair in u_breaks.windows(2) {
318        let half_u = (upair[1] - upair[0]) * 0.5;
319        let middle_u = (upair[1] + upair[0]) * 0.5;
320        for vpair in v_breaks.windows(2) {
321            let half_v = (vpair[1] - vpair[0]) * 0.5;
322            let middle_v = (vpair[1] + vpair[0]) * 0.5;
323            for i in 0..GAUSS_X.len() {
324                for j in 0..GAUSS_X.len() {
325                    total += GAUSS_W[i]
326                        * GAUSS_W[j]
327                        * half_u
328                        * half_v
329                        * evaluate_integrand(
330                            face,
331                            middle_u + half_u * GAUSS_X[i],
332                            middle_v + half_v * GAUSS_X[j],
333                            kind,
334                        )?;
335                }
336            }
337        }
338    }
339    Ok(total)
340}
341
342/// Integrate each of `kinds` over a bi-periodic band face, or return None when
343/// `face` is not such a band. The between-rims strip integrates over the full
344/// periodic span × the cross-level strip; the complement is the full-domain
345/// integral MINUS that strip (the two tile the closed cross period, so no
346/// domain-extended evaluation is needed).
347pub(super) fn biperiodic_band_integral(
348    face: &FaceRecord,
349    kinds: &[Integrand],
350) -> Result<Option<Vec<f64>>, String> {
351    let Some(band) = biperiodic_band_range(face)? else {
352        return Ok(None);
353    };
354    let [u0, u1] = face.surface.domain_u()?;
355    let [v0, v1] = face.surface.domain_v()?;
356    let (u_lo, u_hi, v_lo, v_hi) = if band.p_is_u {
357        (u0, u1, band.q_lo, band.q_hi)
358    } else {
359        (band.q_lo, band.q_hi, v0, v1)
360    };
361    let mut out = Vec::with_capacity(kinds.len());
362    for &kind in kinds {
363        let strip = integrate_rectangle(face, u_lo, u_hi, v_lo, v_hi, kind)?;
364        let value = if band.complement {
365            integrate_untrimmed(face, kind)? - strip
366        } else {
367            strip
368        };
369        out.push(value);
370    }
371    Ok(Some(out))
372}
373
374/// Append one coedge's pcurve samples (shifted by `offset`) to `polygon`.
375pub(super) fn sample_coedge(
376    coedge: &crate::topology::CoedgeRecord,
377    offset: [f64; 2],
378    polygon: &mut Vec<[f64; 2]>,
379) -> Result<(), String> {
380    let curve = &coedge.pcurve;
381    let [start, end] = curve.domain()?;
382    let interior = interior_knots(&curve.knots, curve.degree);
383    let rational = curve
384        .control_points
385        .iter()
386        .any(|point| (point.w - curve.control_points[0].w).abs() > 1e-12);
387    if curve.degree == 1 && !rational {
388        for parameter in std::iter::once(start).chain(interior.iter().copied()) {
389            let point = curve.evaluate(parameter)?;
390            polygon.push([point.x + offset[0], point.y + offset[1]]);
391        }
392    } else {
393        let sample_count = (24usize).max((interior.len() + 1) * 16).min(128);
394        for index in 0..sample_count {
395            let point =
396                curve.evaluate(start + (end - start) * index as f64 / sample_count as f64)?;
397            polygon.push([point.x + offset[0], point.y + offset[1]]);
398        }
399    }
400    Ok(())
401}
brep_kernel/props/mass_properties/integration.rs

brep_kernel/props/mass_properties/
integration.rs
