Skip to main content

geop_core_topology/contains/
face.rs

1use geop_core_geometry::{
2    contains::curve::curve_could_contain,
3    intersection::curve_curve_intersect,
4    nurb_curve::{NurbCurve, NurbCurve2D},
5    nurb_surface::NurbSurface3D,
6};
7use geop_core_math::{
8    geop_error::{GeopError, GeopResult},
9    scalars::Scalar,
10    vector::{Vector2, Vector3},
11};
12
13use crate::{CoedgeId, FaceId, Model, boundary::BoundaryType, contains::rng::Rng};
14
15/// Result of classifying a query point against a face's trimmed boundary.
16#[derive(Clone, Copy, Debug, PartialEq, Eq)]
17pub enum PointClassification {
18    /// The query point coincides with a vertex (shared by two coedges).
19    OnVertex,
20    /// The query point lies on a coedge's pcurve, away from its endpoints.
21    OnCoedge,
22    /// The query point is strictly inside the trimmed boundary.
23    Inside,
24    /// The query point is strictly outside the trimmed boundary (or inside
25    /// a hole).
26    Outside,
27}
28
29const MAX_RAY_ATTEMPTS: usize = 64;
30
31/// Classify `(u, v)` against `face_id`'s trimmed boundary (outer loop minus
32/// holes): [`PointClassification::OnVertex`] / [`PointClassification::OnCoedge`]
33/// if the query point itself coincides with a vertex or lies on a coedge's
34/// pcurve, else [`PointClassification::Inside`]/[`PointClassification::Outside`]
35/// via ray casting in parameter space.
36///
37/// The ray direction is drawn from a seeded PRNG (see [`Rng`]) and retried
38/// (up to a bounded number of attempts) until every crossing it finds lands
39/// strictly inside a coedge's pcurve, away from any vertex — vertex grazes
40/// are ambiguous to count (shared by two coedges) so they're avoided rather
41/// than specially classified. Once such a direction is found, the parity
42/// (even/odd) of its crossing count determines inside/outside; this needs no
43/// normal or winding-direction information, so it works regardless of a
44/// face's loop orientation.
45///
46/// `max_nodes` bounds both the `curve_could_contain` BFS subdivision search
47/// (on-vertex/on-coedge tests) and the `curve_curve_intersect` DFS search
48/// (edge-interior hits). `seed` seeds the direction PRNG.
49pub fn face_contains<S: Scalar>(
50    model: &Model<S>,
51    face_id: FaceId,
52    u: S,
53    v: S,
54    max_nodes: usize,
55    epsilon: S,
56    seed: u64,
57) -> GeopResult<PointClassification> {
58    let face = &model.faces[&face_id];
59    let coedges: Vec<CoedgeId> = model.iterate_face_coedges(face_id).collect();
60    loops_contain(
61        model,
62        &face.surface,
63        &coedges,
64        Vector2::from_array([u, v]),
65        max_nodes,
66        epsilon,
67        seed,
68    )
69}
70
71/// The same classification as [`face_contains`], but against an explicit set
72/// of loops rather than all of a face's.
73///
74/// Exists because "inside this face" and "inside this face's *outer* loop"
75/// are different questions, and validation needs the second: a hole has to
76/// lie within the outer boundary, and asking [`face_contains`] would only
77/// ever answer `OnCoedge` for a point taken from the hole itself. Splitting
78/// the ray casting out here keeps one implementation of it rather than a
79/// second copy that could drift.
80pub fn loops_contain<S: Scalar>(
81    model: &Model<S>,
82    surface: &NurbSurface3D<S>,
83    coedges: &[CoedgeId],
84    query: Vector2<S>,
85    max_nodes: usize,
86    epsilon: S,
87    seed: u64,
88) -> GeopResult<PointClassification> {
89    // Is the query point itself a vertex, or on some coedge's pcurve?
90    // Checked as two full passes (all vertices, then all curves) so
91    // `OnVertex` always takes priority over `OnCoedge` regardless of
92    // coedge iteration order.
93    for &coedge_id in coedges {
94        let pcurve = &model.coedges[&coedge_id].pcurve;
95        let vertex_pt = pcurve.evaluate(pcurve.domain().0)?;
96        if vertex_pt.could_be_equal(&query) {
97            return Ok(PointClassification::OnVertex);
98        }
99    }
100    for &coedge_id in coedges {
101        let pcurve = &model.coedges[&coedge_id].pcurve;
102        if curve_could_contain(pcurve, &query, max_nodes, epsilon)?.is_some() {
103            return Ok(PointClassification::OnCoedge);
104        }
105    }
106
107    let (u_lo, u_hi) = surface.domain_u();
108    let (v_lo, v_hi) = surface.domain_v();
109    let du = u_hi.sub(u_lo);
110    let dv = v_hi.sub(v_lo);
111    let diag = du.mul(du).add(dv.mul(dv)).sqrt()?;
112    let ray_length = diag.mul(S::from_f64(3.0)).add(S::ONE);
113    let t_epsilon = epsilon.div(ray_length)?;
114
115    let mut rng = Rng::new(seed);
116    // Why the most recent direction was given up on — reported if every one
117    // is, since "no clear direction" alone says nothing about the cause.
118    let mut last_rejection = String::new();
119    'attempt: for _ in 0..MAX_RAY_ATTEMPTS {
120        let dir = rng.next_direction2::<S>();
121        let far = query.add(&dir.prod_scalar(ray_length));
122        let ray: NurbCurve2D<S> = NurbCurve::try_new(
123            1,
124            vec![
125                Vector3::from_array([query[0], query[1], S::ONE]),
126                Vector3::from_array([far[0], far[1], S::ONE]),
127            ],
128            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
129        )?;
130
131        for &coedge_id in coedges {
132            let pcurve = &model.coedges[&coedge_id].pcurve;
133            let vertex_pt = pcurve.evaluate(pcurve.domain().0)?;
134            if curve_could_contain(&ray, &vertex_pt, max_nodes, epsilon)?.is_some() {
135                last_rejection = format!(
136                    "ray {ray:?} could pass through coedge {coedge_id}'s start {vertex_pt:?}"
137                );
138                continue 'attempt;
139            }
140        }
141
142        let mut count = 0usize;
143        for &coedge_id in coedges {
144            let pcurve = &model.coedges[&coedge_id].pcurve;
145            let (d0, d1) = pcurve.domain();
146            // An `Err` here means the search exhausted its node budget
147            // before converging — the same ambiguous signal as a vertex/edge
148            // graze (see `curve_curve_intersect`'s own doc comment: budget
149            // exhaustion is itself evidence of a near-tangential or
150            // coincident ray, not a reliable crossing count), so it's
151            // handled the same way: retry with a fresh direction rather than
152            // propagating a hard failure.
153            let hits = match curve_curve_intersect(&ray, pcurve, max_nodes, max_nodes, epsilon) {
154                Ok(hits) => hits.into_vec(),
155                Err(e) => {
156                    last_rejection =
157                        format!("ray {ray:?} x coedge {coedge_id} pcurve {pcurve:?}: {e:?}");
158                    continue 'attempt;
159                }
160            };
161            for (t, mid) in hits {
162                if !t.definitely_greater(t_epsilon) {
163                    continue;
164                }
165                // `curve_curve_intersect` honestly returns the whole
166                // surviving span of its converged leaf, not an arbitrarily
167                // narrowed midpoint — sharpen before comparing against the
168                // pcurve's own endpoints.
169                let mid = mid.midpoint();
170                if !mid.sub(d0).abs().definitely_greater(epsilon)
171                    || !mid.sub(d1).abs().definitely_greater(epsilon)
172                {
173                    // Grazes a vertex despite the check above (numerical
174                    // slop right at the boundary) — retry with a fresh
175                    // direction rather than risk mis-counting it.
176                    last_rejection = format!(
177                        "ray {ray:?} grazes coedge {coedge_id}'s end at t={mid:?} (domain {d0:?}..{d1:?})"
178                    );
179                    continue 'attempt;
180                }
181                count += 1;
182            }
183        }
184        return Ok(if count % 2 == 1 {
185            PointClassification::Inside
186        } else {
187            PointClassification::Outside
188        });
189    }
190    Err(GeopError::new(format!(
191        "loops_contain: could not find a ray direction clear of every vertex after many attempts; \
192         the last one was rejected because {last_rejection}"
193    )))
194}
195
196/// How many interior points [`face_interior_point_where`] offers from one
197/// boundary base point before moving to the next: enough for two genuinely
198/// different points, few enough that rejecting everything stays cheap.
199const POINTS_PER_BASE: usize = 2;
200
201/// How many times [`face_interior_point`] may halve its step before giving
202/// up. Bounds effort only: each halving is another attempt to land inside the
203/// trim, and exhausting them is reported as an error rather than accepted.
204const MAX_HALVINGS: usize = 40;
205
206/// A `(u, v)` strictly inside `face_id`'s trimmed region.
207///
208/// Found the way the operation itself suggests: start on the boundary and
209/// step inward. The step is taken along the inward normal of the outer loop
210/// at a boundary point — inward in `(u, v)`, obtained by rotating the loop's
211/// own tangent — and halved whenever the point it lands on is not
212/// `Inside`. Halving converges on any non-degenerate face, since a
213/// sufficiently short inward step from a boundary point is always interior,
214/// and it needs no guess about the face's size.
215///
216/// The domain midpoint is not usable for this: a trimmed face need not
217/// contain it, and for a face carved out by a boolean's remesh it very often
218/// does not.
219///
220/// `max_nodes`/`epsilon`/`seed` are passed straight to [`face_contains`].
221pub fn face_interior_point<S: Scalar>(
222    model: &Model<S>,
223    face_id: FaceId,
224    max_nodes: usize,
225    epsilon: S,
226    seed: u64,
227) -> GeopResult<(S, S)> {
228    let found =
229        face_interior_point_where(model, face_id, max_nodes, epsilon, seed, |_, _| Ok(true))?;
230    Ok(found.expect("the first interior point found is always accepted"))
231}
232
233/// Like [`face_interior_point`], but for a caller that needs a point with
234/// some further property: up to one interior point per outer coedge (stepped
235/// in from it, in loop order) is handed to `accept`, and the first it accepts
236/// is returned. `Ok(None)` means interior points were found but none was
237/// accepted; an error, as for `face_interior_point`, that none was found.
238pub fn face_interior_point_where<S: Scalar>(
239    model: &Model<S>,
240    face_id: FaceId,
241    max_nodes: usize,
242    epsilon: S,
243    seed: u64,
244    mut accept: impl FnMut(S, S) -> GeopResult<bool>,
245) -> GeopResult<Option<(S, S)>> {
246    let face = &model.faces[&face_id];
247    let BoundaryType::Loop(anchor) = face.outer else {
248        return Err(GeopError::new(format!(
249            "face_interior_point: face {face_id} is bounded by a bare vertex, so it has no interior to sample"
250        )));
251    };
252
253    let (u_lo, u_hi) = face.surface.domain_u();
254    let (v_lo, v_hi) = face.surface.domain_v();
255    let du = u_hi.sub(u_lo);
256    let dv = v_hi.sub(v_lo);
257    let diagonal = if dv.definitely_greater(du) { dv } else { du };
258
259    // Every coedge of the outer loop is a candidate base point, not just the
260    // anchor's. One base point is not enough in practice: a loop can pass
261    // through a degenerate stretch (a revolve pole, a sliver left by a face
262    // split) where the tangent is unusable or where the face is locally
263    // thinner than `face_contains`' own tolerance band, and there the search
264    // fails however finely it steps — while a different side of the very same
265    // face offers an easy interior point.
266    let coedges: Vec<CoedgeId> = model
267        .iterate_loop_coedges(anchor)
268        .take(model.coedges.len() + 1)
269        .collect();
270
271    let mut found_any = false;
272    'base: for &coedge_id in &coedges {
273        let pcurve = &model.get_coedge(coedge_id)?.pcurve;
274        let (t0, t1) = pcurve.domain();
275        let t = t0.add(t1).div(S::TWO)?.sharpen();
276        let Ok(base) = pcurve.evaluate(t) else {
277            continue;
278        };
279        let Ok(tangent) = pcurve.tangent(t).and_then(|d| d.normalize()) else {
280            continue;
281        };
282
283        // Rotate the tangent a quarter turn in `(u, v)`. Which of the two
284        // perpendiculars points *into* the face depends on the loop's
285        // winding, so both are tried and whichever lands inside wins —
286        // cheaper and more robust than deriving the winding.
287        let normals = [
288            Vector2::from_array([tangent[1].neg(), tangent[0]]),
289            Vector2::from_array([tangent[1], tangent[0].neg()]),
290        ];
291
292        let mut step = diagonal.div(S::TWO)?;
293        let mut offered = 0;
294        for _ in 0..MAX_HALVINGS {
295            for inward in normals {
296                // Any point inside the face will do — a free choice — so the
297                // candidate is sharp: it doesn't inherit the width of the
298                // pcurve it was stepped from, and every later computation
299                // from it (a classification ray, say) starts from a point.
300                let u = base[0].add(inward[0].mul(step)).sharpen();
301                let v = base[1].add(inward[1].mul(step)).sharpen();
302                if u.definitely_less(u_lo)
303                    || u.definitely_greater(u_hi)
304                    || v.definitely_less(v_lo)
305                    || v.definitely_greater(v_hi)
306                {
307                    continue;
308                }
309                // Accepted only if the whole box of radius `epsilon` around it
310                // is inside, not just the point: the point is a free choice,
311                // and every later use of it (a classification ray, a surface
312                // evaluation compared against other solids) works at that
313                // resolution. A point merely strictly inside can sit 5e-16
314                // off the boundary — one inward step exactly the face's width
315                // lands there — and at `epsilon` it *is* on the boundary.
316                let neighbourhood = |t: S| t.sub(epsilon).union(t.add(epsilon));
317                if matches!(
318                    face_contains(
319                        model,
320                        face_id,
321                        neighbourhood(u),
322                        neighbourhood(v),
323                        max_nodes,
324                        epsilon,
325                        seed
326                    )?,
327                    PointClassification::Inside
328                ) {
329                    found_any = true;
330                    if accept(u, v)? {
331                        return Ok(Some((u, v)));
332                    }
333                    // Rejected: the next point comes from half this step, so
334                    // it differs from this one — stepping in from each side
335                    // of a symmetric face lands every base on the same
336                    // centre, so bases alone don't give distinct points. Two
337                    // points per base keep a caller that rejects everything
338                    // (every point it tries turns out alike) cheap.
339                    offered += 1;
340                    if offered == POINTS_PER_BASE {
341                        continue 'base;
342                    }
343                    // Leave the direction loop, reaching the halving below.
344                    break;
345                }
346            }
347            // Halving walks toward the boundary, and `face_contains` reports
348            // `OnCoedge` for anything within its own tolerance of it — so
349            // below that tolerance every step is `OnCoedge` and no amount of
350            // further halving can succeed. Stop and let the next base point
351            // try instead of burning the remaining budget here.
352            if !step.definitely_greater(epsilon) {
353                break;
354            }
355            step = step.div(S::TWO)?;
356        }
357    }
358
359    if found_any {
360        return Ok(None);
361    }
362
363    // Report the outer loop's own `(u, v)` extent. A loop that encloses no
364    // area has nothing inside it and the failure is correct — a degenerate
365    // face, to be fixed wherever it was created. A loop with real extent
366    // means the search failed on a face that does have an interior, which is
367    // this function's bug. The two need opposite fixes and are otherwise
368    // indistinguishable from the message.
369    let mut u_extent = None;
370    let mut v_extent = None;
371    for &coedge_id in &coedges {
372        let Ok(coedge) = model.get_coedge(coedge_id) else {
373            continue;
374        };
375        let (t0, t1) = coedge.pcurve.domain();
376        for i in 0..=4 {
377            let Ok(frac) = S::from_ratio(i, 4) else {
378                continue;
379            };
380            let Ok(uv) = coedge.pcurve.evaluate(t0.add(t1.sub(t0).mul(frac))) else {
381                continue;
382            };
383            u_extent = Some(match u_extent {
384                None => uv[0],
385                Some(e) => S::union(e, uv[0]),
386            });
387            v_extent = Some(match v_extent {
388                None => uv[1],
389                Some(e) => S::union(e, uv[1]),
390            });
391        }
392    }
393    Err(GeopError::new(format!(
394        "face_interior_point: no point strictly inside face {face_id} was found, stepping inward from the midpoint of each of its {} outer coedges; that loop spans u={u_extent:?}, v={v_extent:?} (a loop spanning nothing encloses no area, so the face is degenerate)",
395        coedges.len()
396    )))
397}
398
399#[cfg(test)]
400mod interior_point_tests {
401    use super::face_interior_point;
402    use crate::{Model, test_fixtures::test_cube_solid};
403    use geop_core_math::scalars::{ScalInF64, Scalar};
404
405    const MAX: usize = 20000;
406    const SEED: u64 = 99;
407
408    fn eps() -> ScalInF64 {
409        <ScalInF64 as Scalar>::from_f64(1e-4)
410    }
411
412    /// Every face of a plain cube must yield an interior point.
413    #[test]
414    fn cube_faces_all_have_interior_points() {
415        let mut model = Model::<ScalInF64>::new();
416        let solid = test_cube_solid(&mut model);
417        for face_id in model.solid_faces(solid).unwrap() {
418            face_interior_point(&model, face_id, MAX, eps(), SEED)
419                .unwrap_or_else(|e| panic!("face {face_id}: {e}"));
420        }
421    }
422
423    // `sphere_faces_all_have_interior_points` (exercising seam curves and
424    // degenerate poles, unlike the plain cube above) needs a real
425    // `sphere_solid` — that lives in `geop-ops-extrude-revolve`, which
426    // depends on this crate, so it can't be reached from here without
427    // Cargo compiling this crate twice (see `test_fixtures`'s doc
428    // comment). Covered instead by `geop-ops-extrude-revolve::sphere`'s
429    // own tests, which build the same solid `face_interior_point` runs on
430    // here.
431}
432
433#[cfg(test)]
434mod tests {
435    use super::{PointClassification, face_contains};
436    use crate::{
437        Coedge, CoedgeGeometry, CoedgeId, Edge, Face, FaceId, Model, Sense, ShellId, Vertex,
438        VertexId, boundary::BoundaryType, model::Curve3,
439    };
440    use geop_core_geometry::{
441        nurb_curve::{NurbCurve, NurbCurve2D},
442        nurb_surface::NurbSurface3D,
443    };
444    use geop_core_math::{
445        for_all_scalars,
446        scalars::Scalar,
447        vector::{Vector3, Vector4},
448    };
449
450    const MAX: usize = 200;
451    const EPS: f64 = 1e-3;
452    const SEED: u64 = 12345;
453
454    fn p2<S: Scalar>(x: f64, y: f64) -> Vector3<S> {
455        Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ONE])
456    }
457
458    fn line2<S: Scalar>(a: (f64, f64), b: (f64, f64)) -> NurbCurve2D<S> {
459        NurbCurve::try_new(
460            1,
461            vec![p2(a.0, a.1), p2(b.0, b.1)],
462            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
463        )
464        .unwrap()
465    }
466
467    /// A polygon face on the unit-square `[0,1]^2` parameter surface, built
468    /// from `points` (CCW, `(u, v) == (x, y)`).
469    fn polygon_face<S: Scalar>(model: &mut Model<S>, points: &[(f64, f64)]) -> FaceId {
470        let p =
471            |x: f64, y: f64| Vector4::from_array([S::from_f64(x), S::from_f64(y), S::ZERO, S::ONE]);
472        let surface = NurbSurface3D::try_new(
473            1,
474            1,
475            vec![p(0.0, 0.0), p(0.0, 1.0), p(1.0, 0.0), p(1.0, 1.0)],
476            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
477            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
478        )
479        .unwrap();
480
481        let face_id = model.insert_face(Face {
482            surface,
483            outer: BoundaryType::Vertex(VertexId(0)),
484            holes: Vec::new(),
485            shell: ShellId(999),
486        });
487
488        let n = points.len();
489        let verts: Vec<VertexId> = points
490            .iter()
491            .map(|&(x, y)| {
492                model.insert_vertex(Vertex {
493                    point: Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ZERO]),
494                })
495            })
496            .collect();
497        let edges = (0..n)
498            .map(|i| {
499                model.insert_edge(Edge {
500                    curve: Curve3::try_new(
501                        1,
502                        vec![
503                            Vector4::from_array([
504                                S::from_f64(points[i].0),
505                                S::from_f64(points[i].1),
506                                S::ZERO,
507                                S::ONE,
508                            ]),
509                            Vector4::from_array([
510                                S::from_f64(points[(i + 1) % n].0),
511                                S::from_f64(points[(i + 1) % n].1),
512                                S::ZERO,
513                                S::ONE,
514                            ]),
515                        ],
516                        vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
517                    )
518                    .unwrap(),
519                    start_vertex: verts[i],
520                    end_vertex: verts[(i + 1) % n],
521                })
522            })
523            .collect::<Vec<_>>();
524        let coedges: Vec<CoedgeId> = (0..n)
525            .map(|i| {
526                model.insert_coedge(Coedge {
527                    geometry: CoedgeGeometry::Edge(edges[i]),
528                    sense: Sense::Forward,
529                    pcurve: line2(points[i], points[(i + 1) % n]),
530                    next: CoedgeId(0),
531                    prev: CoedgeId(0),
532                    face: face_id,
533                })
534            })
535            .collect();
536        for i in 0..n {
537            model.coedges.get_mut(&coedges[i]).unwrap().next = coedges[(i + 1) % n];
538            model.coedges.get_mut(&coedges[i]).unwrap().prev = coedges[(i + n - 1) % n];
539        }
540        model.faces.get_mut(&face_id).unwrap().outer = BoundaryType::Loop(coedges[0]);
541
542        face_id
543    }
544
545    /// Diamond with corners at (1,.5) right, (.5,1) top, (0,.5) left,
546    /// (.5,0) bottom, traversed CCW.
547    fn diamond_face<S: Scalar>(model: &mut Model<S>) -> FaceId {
548        polygon_face(model, &[(1., 0.5), (0.5, 1.), (0., 0.5), (0.5, 0.)])
549    }
550
551    fn check_diamond_interior_point_is_contained<S: Scalar>() {
552        let mut model = Model::<S>::new();
553        let face_id = diamond_face(&mut model);
554        assert_eq!(
555            face_contains(
556                &model,
557                face_id,
558                S::from_f64(0.5),
559                S::from_f64(0.3),
560                MAX,
561                S::from_f64(EPS),
562                SEED
563            )
564            .unwrap(),
565            PointClassification::Inside
566        );
567    }
568    #[test]
569    fn diamond_interior_point_is_contained() {
570        for_all_scalars!(check_diamond_interior_point_is_contained);
571    }
572
573    fn check_diamond_exterior_point_is_not_contained<S: Scalar>() {
574        let mut model = Model::<S>::new();
575        let face_id = diamond_face(&mut model);
576        assert_eq!(
577            face_contains(
578                &model,
579                face_id,
580                S::from_f64(0.1),
581                S::from_f64(0.3),
582                MAX,
583                S::from_f64(EPS),
584                SEED
585            )
586            .unwrap(),
587            PointClassification::Outside
588        );
589    }
590    #[test]
591    fn diamond_exterior_point_is_not_contained() {
592        for_all_scalars!(check_diamond_exterior_point_is_not_contained);
593    }
594
595    fn check_diamond_center_hits_convex_vertex_from_inside<S: Scalar>() {
596        let mut model = Model::<S>::new();
597        let face_id = diamond_face(&mut model);
598        assert_eq!(
599            face_contains(
600                &model,
601                face_id,
602                S::from_f64(0.5),
603                S::from_f64(0.5),
604                MAX,
605                S::from_f64(EPS),
606                SEED
607            )
608            .unwrap(),
609            PointClassification::Inside
610        );
611    }
612    #[test]
613    fn diamond_center_hits_convex_vertex_from_inside() {
614        for_all_scalars!(check_diamond_center_hits_convex_vertex_from_inside);
615    }
616
617    fn check_diamond_vertex_query_is_on_vertex<S: Scalar>() {
618        let mut model = Model::<S>::new();
619        let face_id = diamond_face(&mut model);
620        assert_eq!(
621            face_contains(
622                &model,
623                face_id,
624                S::ONE,
625                S::from_f64(0.5),
626                MAX,
627                S::from_f64(EPS),
628                SEED
629            )
630            .unwrap(),
631            PointClassification::OnVertex
632        );
633    }
634    #[test]
635    fn diamond_vertex_query_is_on_vertex() {
636        for_all_scalars!(check_diamond_vertex_query_is_on_vertex);
637    }
638
639    fn check_diamond_edge_query_is_on_coedge<S: Scalar>() {
640        let mut model = Model::<S>::new();
641        let face_id = diamond_face(&mut model);
642        assert_eq!(
643            face_contains(
644                &model,
645                face_id,
646                S::from_f64(0.75),
647                S::from_f64(0.75),
648                MAX,
649                S::from_f64(EPS),
650                SEED
651            )
652            .unwrap(),
653            PointClassification::OnCoedge
654        );
655    }
656    #[test]
657    fn diamond_edge_query_is_on_coedge() {
658        for_all_scalars!(check_diamond_edge_query_is_on_coedge);
659    }
660}