--- a/src/fe_model/p2_local_model.rs Thu Sep 03 10:50:27 2026 -0500 +++ b/src/fe_model/p2_local_model.rs Thu Sep 03 10:59:13 2026 -0500 @@ -8,7 +8,7 @@ use crate::linsolve::*; use crate::loc::Loc; use crate::sets::Cube; -use crate::sets::{NPolygon, Set, SpannedHalfspace}; +use crate::sets::Set; use crate::types::*; use numeric_literals::replace_float_literals; @@ -102,13 +102,13 @@ #[inline] fn p2powers(&self) -> Self::Output { let &Loc([x0, x1]) = self; - [x0 * x0, x0 * x1, x1 * x1].into() + [x0 * x0, 2.0 * x0 * x1, x1 * x1].into() } #[inline] fn p2powers_full(&self) -> Self::Full { let &Loc([x0, x1]) = self; - [1.0, x0, x1, x0 * x0, x0 * x1, x1 * x1].into() + [1.0, x0, x1, x0 * x0, 2.0 * x0 * x1, x1 * x1].into() } #[inline] @@ -328,13 +328,14 @@ let &P2LocalModel { a0, a1: Loc([a1, a2]), a2: Loc([a11, a12, a22]), .. } = self; let &Loc([x00, x01]) = x0; let d @ Loc([d0, d1]) = x1 - x0; - let b0 = a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + a12 * x00 * x01 + a22 * x01 * x01; + let b0 = + a0 + a1 * x00 + a2 * x01 + a11 * x00 * x00 + 2.0 * a12 * x00 * x01 + a22 * x01 * x01; let b1 = a1 * d0 + a2 * d1 + 2.0 * a11 * d0 * x00 - + a12 * (d0 * x01 + d1 * x00) + + 2.0 * a12 * (d0 * x01 + d1 * x00) + 2.0 * a22 * d1 * x01; - let b11 = a11 * d0 * d0 + a12 * d0 * d1 + a22 * d1 * d1; + let b11 = a11 * d0 * d0 + 2.0 * a12 * d0 * d1 + a22 * d1 * d1; let edge_1d_model = P2LocalModel { a0: b0, a1: Loc([b1]), @@ -364,9 +365,9 @@ let r = 2.0 * (a11 * a22 - a12 * a12); if r > 0.0 { // An interior solution (x[1], x[2]) has to satisfy - // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[1] + 2a₁₂*x[1]+a₂=0 + // 2a₁₁*x[1] + 2a₁₂*x[2]+a₁ =0 and 2a₂₂*x[2] + 2a₁₂*x[1]+a₂=0 // This gives - let x = [(a22 * a1 - a12 * a2) / r, (a12 * a1 - a11 * a2) / r].into(); + let x = [(a22 * a1 - a12 * a2) / r, (a11 * a2 - a12 * a1) / r].into(); if el.contains(&x) { return (x, self.value(&x)); }