diff -r 30af49ae2135 -r ce37ff3ce507 src/fe_model/p2_local_model.rs
--- a/src/fe_model/p2_local_model.rs	Sun Aug 23 22:12:08 2026 -0500
+++ b/src/fe_model/p2_local_model.rs	Thu Sep 03 10:14:13 2026 -0500
@@ -39,16 +39,23 @@
     }
 }
 
+#[replace_float_literals(F::cast_from(literal))]
 impl<'a, F: Float> Set<Loc<2, F>> for PlanarSimplex<F> {
     #[inline]
     fn contains<I: Instance<Loc<2, F>>>(&self, z: I) -> bool {
-        let &[x0, x1, x2] = &self.0;
-        NPolygon([
-            [x0, x1].spanned_halfspace(),
-            [x1, x2].spanned_halfspace(),
-            [x2, x0].spanned_halfspace(),
-        ])
-        .contains(z)
+        z.eval_ref(|&Loc([x, y])| {
+            let [Loc([x1, y1]), Loc([x2, y2]), Loc([x3, y3])] = self.0;
+
+            let areax2 = (x1 - x3) * (y2 - y3) + (x2 - x3) * (y1 - y3);
+
+            //  Unscaled barycentric coordinates
+            let s_unscaled = (y2 - y3) * (x - x3) + (x3 - x2) * (y - y3);
+            let t_unscaled = (y3 - y1) * (x - x3) + (x1 - x3) * (y - y3);
+
+            (0.0 <= s_unscaled/*&& s <= areax2*/)
+                && (0.0 <= t_unscaled/*&& t <= areax2*/)
+                && s_unscaled + t_unscaled <= areax2
+        })
     }
 }