Skip to main content

brep_kernel/meshing/
tessellation.rs

1use crate::mass_properties::trim_polygons;
2use crate::topology::{BrepSolid, FaceRecord};
3use crate::{KnotVector, Mesh, Vec3};
4
5#[derive(Clone, Copy, Debug)]
6pub struct TessellationOptions {
7    pub slabs_per_span_u: usize,
8    pub steps_per_span_v: usize,
9}
10
11impl Default for TessellationOptions {
12    fn default() -> Self {
13        Self {
14            slabs_per_span_u: 8,
15            steps_per_span_v: 8,
16        }
17    }
18}
19
20fn interior_knots(knots: &[f64], degree: usize) -> Vec<f64> {
21    let start = knots[degree];
22    let end = knots[knots.len() - 1 - degree];
23    let mut result = Vec::new();
24    for &knot in knots {
25        if knot <= start + 1e-12 || knot >= end - 1e-12 {
26            continue;
27        }
28        if result
29            .last()
30            .is_none_or(|previous: &f64| (knot - *previous).abs() > 1e-12)
31        {
32            result.push(knot);
33        }
34    }
35    result
36}
37
38fn crossings_at(polygons: &[Vec<[f64; 2]>], u: f64) -> Vec<f64> {
39    let mut crossings = Vec::new();
40    for polygon in polygons {
41        for index in 0..polygon.len() {
42            let a = polygon[index];
43            let b = polygon[(index + 1) % polygon.len()];
44            if (a[0] > u) != (b[0] > u) {
45                crossings.push(a[1] + (u - a[0]) / (b[0] - a[0]) * (b[1] - a[1]));
46            }
47        }
48    }
49    crossings.sort_by(f64::total_cmp);
50    crossings
51}
52
53fn emit_vertex(
54    mesh: &mut Mesh,
55    face: &FaceRecord,
56    u: f64,
57    v: f64,
58    domains: [f64; 4],
59) -> Result<u32, String> {
60    let (point, su, sv) = face.surface.deriv1(u, v)?;
61    let mut normal = su.cross(sv);
62    if normal.length() <= 1e-12 {
63        let [u0, u1, v0, v1] = domains;
64        let epsilon = 1e-4;
65        let nudged_u = u.clamp(u0 + epsilon * (u1 - u0), u1 - epsilon * (u1 - u0));
66        let nudged_v = v.clamp(v0 + epsilon * (v1 - v0), v1 - epsilon * (v1 - v0));
67        let (_, nsu, nsv) = face.surface.deriv1(nudged_u, nudged_v)?;
68        normal = nsu.cross(nsv);
69        if normal.length() <= 1e-12 {
70            normal = Vec3::new(0.0, 0.0, 1.0);
71        }
72    }
73    normal = normal.normalized()?;
74    if !face.same_sense {
75        normal = normal.scale(-1.0);
76    }
77    mesh.positions.extend([point.x, point.y, point.z]);
78    mesh.normals.extend([normal.x, normal.y, normal.z]);
79    Ok((mesh.positions.len() / 3 - 1) as u32)
80}
81
82fn push_triangle(mesh: &mut Mesh, a: u32, b: u32, c: u32, ccw: bool, face_id: u32) {
83    if ccw {
84        mesh.indices.extend([a, b, c]);
85    } else {
86        mesh.indices.extend([a, c, b]);
87    }
88    mesh.face_ids.push(face_id);
89}
90
91pub fn tessellate_face(
92    face: &FaceRecord,
93    options: TessellationOptions,
94    face_id: u32,
95) -> Result<Mesh, String> {
96    let ku = KnotVector::new(face.surface.knots_u.clone(), face.surface.degree_u)?;
97    let kv = KnotVector::new(face.surface.knots_v.clone(), face.surface.degree_v)?;
98    let [u0, u1] = ku.domain();
99    let [v0, v1] = kv.domain();
100    let polygons = trim_polygons(face)?;
101    let slabs = options.slabs_per_span_u.max(1);
102    let steps = options.steps_per_span_v.max(1);
103
104    let mut u_knots = vec![u0];
105    u_knots.extend(interior_knots(&face.surface.knots_u, face.surface.degree_u));
106    u_knots.push(u1);
107    let mut breaks = Vec::new();
108    for pair in u_knots.windows(2) {
109        for index in 0..=slabs {
110            breaks.push(pair[0] + (pair[1] - pair[0]) * index as f64 / slabs as f64);
111        }
112    }
113    for polygon in &polygons {
114        breaks.extend(polygon.iter().map(|point| point[0].clamp(u0, u1)));
115    }
116    breaks.sort_by(f64::total_cmp);
117    let minimum_gap = (u1 - u0) * 1e-9;
118    breaks.dedup_by(|a, b| (*a - *b).abs() <= minimum_gap);
119
120    let v_span_length =
121        (v1 - v0) / (interior_knots(&face.surface.knots_v, face.surface.degree_v).len() + 1) as f64;
122    let v_step = v_span_length / steps as f64;
123    let mut mesh = Mesh::default();
124    for pair in breaks.windows(2) {
125        let ua = pair[0];
126        let ub = pair[1];
127        if ub - ua <= minimum_gap {
128            continue;
129        }
130        let midpoint_crossings = crossings_at(&polygons, (ua + ub) * 0.5);
131        if midpoint_crossings.len() < 2 {
132            continue;
133        }
134        let mut left_crossings = crossings_at(&polygons, ua + (ub - ua) * 1e-7);
135        let mut right_crossings = crossings_at(&polygons, ub - (ub - ua) * 1e-7);
136        if left_crossings.len() != midpoint_crossings.len()
137            || right_crossings.len() != midpoint_crossings.len()
138        {
139            left_crossings.clone_from(&midpoint_crossings);
140            right_crossings.clone_from(&midpoint_crossings);
141        }
142        for interval in (0..midpoint_crossings.len() - 1).step_by(2) {
143            let lower_left = left_crossings[interval].clamp(v0, v1);
144            let upper_left = left_crossings[interval + 1].clamp(v0, v1);
145            let lower_right = right_crossings[interval].clamp(v0, v1);
146            let upper_right = right_crossings[interval + 1].clamp(v0, v1);
147            let span = (upper_left - lower_left).max(upper_right - lower_right);
148            if span <= 1e-12 {
149                continue;
150            }
151            let row_steps = ((span / v_step).ceil() as usize).clamp(1, 64);
152            let mut left = Vec::with_capacity(row_steps + 1);
153            let mut right = Vec::with_capacity(row_steps + 1);
154            for index in 0..=row_steps {
155                let fraction = index as f64 / row_steps as f64;
156                left.push(emit_vertex(
157                    &mut mesh,
158                    face,
159                    ua,
160                    lower_left + (upper_left - lower_left) * fraction,
161                    [u0, u1, v0, v1],
162                )?);
163                right.push(emit_vertex(
164                    &mut mesh,
165                    face,
166                    ub,
167                    lower_right + (upper_right - lower_right) * fraction,
168                    [u0, u1, v0, v1],
169                )?);
170            }
171            for index in 0..row_steps {
172                push_triangle(
173                    &mut mesh,
174                    left[index],
175                    right[index],
176                    right[index + 1],
177                    face.same_sense,
178                    face_id,
179                );
180                push_triangle(
181                    &mut mesh,
182                    left[index],
183                    right[index + 1],
184                    left[index + 1],
185                    face.same_sense,
186                    face_id,
187                );
188            }
189        }
190    }
191    mesh.validate()?;
192    Ok(mesh)
193}
194
195pub fn tessellate_brep(solid: &BrepSolid, options: TessellationOptions) -> Result<Mesh, String> {
196    let mut mesh = Mesh::default();
197    let mut sequential_face_id = 0u32;
198    for shell in &solid.shells {
199        for face in &shell.faces {
200            let face_mesh = tessellate_face(face, options, sequential_face_id)?;
201            let base = (mesh.positions.len() / 3) as u32;
202            mesh.positions.extend(face_mesh.positions);
203            mesh.normals.extend(face_mesh.normals);
204            mesh.indices
205                .extend(face_mesh.indices.into_iter().map(|index| index + base));
206            mesh.face_ids.extend(face_mesh.face_ids);
207            sequential_face_id += 1;
208        }
209    }
210    mesh.validate()?;
211    Ok(mesh)
212}
213
214#[cfg(test)]
215mod tests {
216    use super::*;
217    use crate::{make_box_brep, make_sphere_brep};
218
219    #[test]
220    fn exact_box_tessellation_has_outward_volume() {
221        let solid = make_box_brep(Vec3::default(), 2.0, 3.0, 4.0).unwrap();
222        let mesh = tessellate_brep(&solid, TessellationOptions::default()).unwrap();
223        assert!((mesh.signed_volume() - 24.0).abs() < 1e-8);
224        assert_eq!(mesh.face_ids.len(), mesh.indices.len() / 3);
225    }
226
227    #[test]
228    fn sphere_tessellation_approximates_analytic_volume() {
229        let solid = make_sphere_brep(Vec3::default(), 3.0, Vec3::new(0.0, 1.0, 0.0)).unwrap();
230        let mesh = tessellate_brep(
231            &solid,
232            TessellationOptions {
233                slabs_per_span_u: 24,
234                steps_per_span_v: 24,
235            },
236        )
237        .unwrap();
238        let expected = 4.0 * std::f64::consts::PI * 27.0 / 3.0;
239        assert!((mesh.signed_volume() - expected).abs() / expected < 0.01);
240    }
241}