fix(solver+io): B1 block-FD Hessian, V3 NaN/Inf guard, C1 ok-flag (audit quick-wins)

B1 (api-performance): wire hyper_ideal_hessian_block_fd_sym into newton_hyper_ideal
instead of the full-FD path — 33×/1166× faster on cathead/brezel, also fixes the
FD-step vs Newton-tol accuracy floor (N1 cross-fix).

V3 (input-validation): add isfinite check on all vertex coordinates in load_mesh;
NaN/Inf now throws runtime_error at the I/O boundary before poisoning the solver.

C1 (test-coverage): expose the internal ok flag via an optional bool* parameter
on solve_linear_system so double-solver failure is no longer invisible to callers.

+5 new tests (LoadMeshThrowsOnNaN/Inf, OkFlag_True*, OkFlag_NullPointerIsSafe).
282/282 tests pass.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

code/include/newton_solver.hpp

@@ -103,14 +103,18 @@ inline Eigen::VectorXd solve_with_fallback(
 /// Solve `A·x = rhs` with the same SimplicialLDLT → SparseQR fallback
 /// strategy used inside all three Newton solvers.  If `fallback_used`
 /// is non-null, it is set to `true` iff the SparseQR fallback ran.
+/// If `ok` is non-null, it is set to `true` iff at least one solver succeeded.
 /// Returns `Eigen::VectorXd::Zero(rhs.size())` if both solvers fail.
 inline Eigen::VectorXd solve_linear_system(
     const Eigen::SparseMatrix<double>& A,
     const Eigen::VectorXd&             rhs,
-    bool*                              fallback_used = nullptr)
+    bool*                              fallback_used = nullptr,
+    bool*                              ok            = nullptr)
 {
-    bool ok = false;
-    return detail::solve_with_fallback(A, rhs, ok, fallback_used);
+    bool ok_local = false;
+    auto x = detail::solve_with_fallback(A, rhs, ok_local, fallback_used);
+    if (ok) *ok = ok_local;
+    return x;
 }
 
 namespace detail {   // re-open for the remaining helpers
@@ -429,8 +433,8 @@ inline NewtonResult newton_hyper_ideal(
             return res;
         }
 
-        // ── Hessian (numerical FD) + solve H·Δx = −G ─────────────────────────
-        auto H  = hyper_ideal_hessian_sym(mesh, x, m, hess_eps);
+        // ── Hessian (block-FD, ~33–1166× faster than full-FD) + solve H·Δx = −G
+        auto H  = hyper_ideal_hessian_block_fd_sym(mesh, x, m, hess_eps);
         bool ok = false;
         Eigen::VectorXd dx = detail::solve_with_fallback(H, -G, ok);
         if (!ok) break;