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

brep_kernel/props/mass_properties/
integration.rs
