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    // Each coedge's edge, start vertex and pcurve start: enough to recognise
394    // the face in the model, to see a loop that doubles back on itself (an
395    // edge appearing twice, as a spur), and to see a pcurve whose `(u, v)`
396    // disagrees with where its 3-D vertex actually is on the surface.
397    let loop_description: Vec<String> = coedges
398        .iter()
399        .map(|&coedge_id| {
400            let start = model.coedge_start_vertex(coedge_id).map(|v| v.point);
401            let coedge = model.get_coedge(coedge_id);
402            let geometry = coedge.as_ref().ok().map(|c| c.geometry);
403            let start_uv = coedge.and_then(|c| c.pcurve.evaluate(c.pcurve.domain().0));
404            format!("{coedge_id} ({geometry:?}): starts at {start:?}, (u, v) = {start_uv:?}")
405        })
406        .collect();
407    Err(GeopError::new(format!(
408        "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) within the surface's domain u={:?}, v={:?}; the loop: [{}]",
409        coedges.len(),
410        (u_lo, u_hi),
411        (v_lo, v_hi),
412        loop_description.join("; ")
413    )))
414}
415
416#[cfg(test)]
417mod interior_point_tests {
418    use super::face_interior_point;
419    use crate::{Model, test_fixtures::test_cube_solid};
420    use geop_core_math::scalars::{ScalInF64, Scalar};
421
422    const MAX: usize = 20000;
423    const SEED: u64 = 99;
424
425    fn eps() -> ScalInF64 {
426        <ScalInF64 as Scalar>::from_f64(1e-4)
427    }
428
429    /// Every face of a plain cube must yield an interior point.
430    #[test]
431    fn cube_faces_all_have_interior_points() {
432        let mut model = Model::<ScalInF64>::new();
433        let solid = test_cube_solid(&mut model);
434        for face_id in model.solid_faces(solid).unwrap() {
435            face_interior_point(&model, face_id, MAX, eps(), SEED)
436                .unwrap_or_else(|e| panic!("face {face_id}: {e}"));
437        }
438    }
439
440    // `sphere_faces_all_have_interior_points` (exercising seam curves and
441    // degenerate poles, unlike the plain cube above) needs a real
442    // `sphere_solid` — that lives in `geop-ops-extrude-revolve`, which
443    // depends on this crate, so it can't be reached from here without
444    // Cargo compiling this crate twice (see `test_fixtures`'s doc
445    // comment). Covered instead by `geop-ops-extrude-revolve::sphere`'s
446    // own tests, which build the same solid `face_interior_point` runs on
447    // here.
448}
449
450#[cfg(test)]
451mod tests {
452    use super::{PointClassification, face_contains};
453    use crate::{
454        Coedge, CoedgeGeometry, CoedgeId, Edge, Face, FaceId, Model, Sense, ShellId, Vertex,
455        VertexId, boundary::BoundaryType, model::Curve3,
456    };
457    use geop_core_geometry::{
458        nurb_curve::{NurbCurve, NurbCurve2D},
459        nurb_surface::NurbSurface3D,
460    };
461    use geop_core_math::{
462        for_all_scalars,
463        scalars::Scalar,
464        vector::{Vector3, Vector4},
465    };
466
467    const MAX: usize = 200;
468    const EPS: f64 = 1e-3;
469    const SEED: u64 = 12345;
470
471    fn p2<S: Scalar>(x: f64, y: f64) -> Vector3<S> {
472        Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ONE])
473    }
474
475    fn line2<S: Scalar>(a: (f64, f64), b: (f64, f64)) -> NurbCurve2D<S> {
476        NurbCurve::try_new(
477            1,
478            vec![p2(a.0, a.1), p2(b.0, b.1)],
479            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
480        )
481        .unwrap()
482    }
483
484    /// A polygon face on the unit-square `[0,1]^2` parameter surface, built
485    /// from `points` (CCW, `(u, v) == (x, y)`).
486    fn polygon_face<S: Scalar>(model: &mut Model<S>, points: &[(f64, f64)]) -> FaceId {
487        let p =
488            |x: f64, y: f64| Vector4::from_array([S::from_f64(x), S::from_f64(y), S::ZERO, S::ONE]);
489        let surface = NurbSurface3D::try_new(
490            1,
491            1,
492            vec![p(0.0, 0.0), p(0.0, 1.0), p(1.0, 0.0), p(1.0, 1.0)],
493            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
494            vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
495        )
496        .unwrap();
497
498        let face_id = model.insert_face(Face {
499            surface,
500            outer: BoundaryType::Vertex(VertexId(0)),
501            holes: Vec::new(),
502            shell: ShellId(999),
503        });
504
505        let n = points.len();
506        let verts: Vec<VertexId> = points
507            .iter()
508            .map(|&(x, y)| {
509                model.insert_vertex(Vertex {
510                    point: Vector3::from_array([S::from_f64(x), S::from_f64(y), S::ZERO]),
511                })
512            })
513            .collect();
514        let edges = (0..n)
515            .map(|i| {
516                model.insert_edge(Edge {
517                    curve: Curve3::try_new(
518                        1,
519                        vec![
520                            Vector4::from_array([
521                                S::from_f64(points[i].0),
522                                S::from_f64(points[i].1),
523                                S::ZERO,
524                                S::ONE,
525                            ]),
526                            Vector4::from_array([
527                                S::from_f64(points[(i + 1) % n].0),
528                                S::from_f64(points[(i + 1) % n].1),
529                                S::ZERO,
530                                S::ONE,
531                            ]),
532                        ],
533                        vec![S::ZERO, S::ZERO, S::ONE, S::ONE],
534                    )
535                    .unwrap(),
536                    start_vertex: verts[i],
537                    end_vertex: verts[(i + 1) % n],
538                })
539            })
540            .collect::<Vec<_>>();
541        let coedges: Vec<CoedgeId> = (0..n)
542            .map(|i| {
543                model.insert_coedge(Coedge {
544                    geometry: CoedgeGeometry::Edge(edges[i]),
545                    sense: Sense::Forward,
546                    pcurve: line2(points[i], points[(i + 1) % n]),
547                    next: CoedgeId(0),
548                    prev: CoedgeId(0),
549                    face: face_id,
550                })
551            })
552            .collect();
553        for i in 0..n {
554            model.coedges.get_mut(&coedges[i]).unwrap().next = coedges[(i + 1) % n];
555            model.coedges.get_mut(&coedges[i]).unwrap().prev = coedges[(i + n - 1) % n];
556        }
557        model.faces.get_mut(&face_id).unwrap().outer = BoundaryType::Loop(coedges[0]);
558
559        face_id
560    }
561
562    /// Diamond with corners at (1,.5) right, (.5,1) top, (0,.5) left,
563    /// (.5,0) bottom, traversed CCW.
564    fn diamond_face<S: Scalar>(model: &mut Model<S>) -> FaceId {
565        polygon_face(model, &[(1., 0.5), (0.5, 1.), (0., 0.5), (0.5, 0.)])
566    }
567
568    fn check_diamond_interior_point_is_contained<S: Scalar>() {
569        let mut model = Model::<S>::new();
570        let face_id = diamond_face(&mut model);
571        assert_eq!(
572            face_contains(
573                &model,
574                face_id,
575                S::from_f64(0.5),
576                S::from_f64(0.3),
577                MAX,
578                S::from_f64(EPS),
579                SEED
580            )
581            .unwrap(),
582            PointClassification::Inside
583        );
584    }
585    #[test]
586    fn diamond_interior_point_is_contained() {
587        for_all_scalars!(check_diamond_interior_point_is_contained);
588    }
589
590    fn check_diamond_exterior_point_is_not_contained<S: Scalar>() {
591        let mut model = Model::<S>::new();
592        let face_id = diamond_face(&mut model);
593        assert_eq!(
594            face_contains(
595                &model,
596                face_id,
597                S::from_f64(0.1),
598                S::from_f64(0.3),
599                MAX,
600                S::from_f64(EPS),
601                SEED
602            )
603            .unwrap(),
604            PointClassification::Outside
605        );
606    }
607    #[test]
608    fn diamond_exterior_point_is_not_contained() {
609        for_all_scalars!(check_diamond_exterior_point_is_not_contained);
610    }
611
612    fn check_diamond_center_hits_convex_vertex_from_inside<S: Scalar>() {
613        let mut model = Model::<S>::new();
614        let face_id = diamond_face(&mut model);
615        assert_eq!(
616            face_contains(
617                &model,
618                face_id,
619                S::from_f64(0.5),
620                S::from_f64(0.5),
621                MAX,
622                S::from_f64(EPS),
623                SEED
624            )
625            .unwrap(),
626            PointClassification::Inside
627        );
628    }
629    #[test]
630    fn diamond_center_hits_convex_vertex_from_inside() {
631        for_all_scalars!(check_diamond_center_hits_convex_vertex_from_inside);
632    }
633
634    fn check_diamond_vertex_query_is_on_vertex<S: Scalar>() {
635        let mut model = Model::<S>::new();
636        let face_id = diamond_face(&mut model);
637        assert_eq!(
638            face_contains(
639                &model,
640                face_id,
641                S::ONE,
642                S::from_f64(0.5),
643                MAX,
644                S::from_f64(EPS),
645                SEED
646            )
647            .unwrap(),
648            PointClassification::OnVertex
649        );
650    }
651    #[test]
652    fn diamond_vertex_query_is_on_vertex() {
653        for_all_scalars!(check_diamond_vertex_query_is_on_vertex);
654    }
655
656    fn check_diamond_edge_query_is_on_coedge<S: Scalar>() {
657        let mut model = Model::<S>::new();
658        let face_id = diamond_face(&mut model);
659        assert_eq!(
660            face_contains(
661                &model,
662                face_id,
663                S::from_f64(0.75),
664                S::from_f64(0.75),
665                MAX,
666                S::from_f64(EPS),
667                SEED
668            )
669            .unwrap(),
670            PointClassification::OnCoedge
671        );
672    }
673    #[test]
674    fn diamond_edge_query_is_on_coedge() {
675        for_all_scalars!(check_diamond_edge_query_is_on_coedge);
676    }
677}