diff --git a/src/cvc/geometry/cvc-mesher/LBIE/octree.cpp b/src/cvc/geometry/cvc-mesher/LBIE/octree.cpp
index 0704fd16b..f50f6f5e0 100644
--- a/src/cvc/geometry/cvc-mesher/LBIE/octree.cpp
+++ b/src/cvc/geometry/cvc-mesher/LBIE/octree.cpp
@@ -5924,18 +5924,20 @@ void Octree::getGrad(float v0[3], float v1[3], float v2[3], float v3[3], float&
 	else jacob = (float) sqrt(jacob);
 	condition_no = (g11 + g22 + g33) / (3.0f*jacob);
 
+	// t1 is d(jacobian)/d(v0) for the edges v1-v0, v2-v0, v3-v0 (unnormalized). Its second factors
+	// were indexed 1-based: [3] for z, which read past these float[3] arrays, and [2] for y.
 	t0 = 2.0f*(v1[0] + v2[0] + v3[0] - 3.0f*v0[0]);
-	t1 = (v1[1]-v0[1])*(v3[3]-v2[3]) + (v2[1]-v0[1])*(v1[3]-v3[3]) + (v3[1]-v0[1])*(v2[3]-v1[3]);
+	t1 = (v1[1]-v0[1])*(v3[2]-v2[2]) + (v2[1]-v0[1])*(v1[2]-v3[2]) + (v3[1]-v0[1])*(v2[2]-v1[2]);
 	grad[0] = t0/jacob - (g11+g22+g33)*t1/(jacob*jacob);
 	//grad[0] = t1/(g11+g22+g33) - jacob*t0/((g11+g22+g33)*(g11+g22+g33));
 
 	t0 = 2.0f*(v1[1] + v2[1] + v3[1] - 3.0f*v0[1]);
-	t1 = (v1[0]-v0[0])*(v2[3]-v3[3]) + (v2[0]-v0[0])*(v3[3]-v1[3]) + (v3[0]-v0[0])*(v1[3]-v2[3]);
+	t1 = (v1[0]-v0[0])*(v2[2]-v3[2]) + (v2[0]-v0[0])*(v3[2]-v1[2]) + (v3[0]-v0[0])*(v1[2]-v2[2]);
 	grad[1] = t0/jacob - (g11+g22+g33)*t1/(jacob*jacob);
 	//grad[1] = t1/(g11+g22+g33) - jacob*t0/((g11+g22+g33)*(g11+g22+g33));
 
 	t0 = 2.0f*(v1[2] + v2[2] + v3[2] - 3.0f*v0[2]);
-	t1 = (v1[0]-v0[0])*(v3[2]-v2[2]) + (v2[0]-v0[0])*(v1[2]-v3[2]) + (v3[0]-v0[0])*(v2[2]-v1[2]);
+	t1 = (v1[0]-v0[0])*(v3[1]-v2[1]) + (v2[0]-v0[0])*(v1[1]-v3[1]) + (v3[0]-v0[0])*(v2[1]-v1[1]);
 	grad[2] = t0/jacob - (g11+g22+g33)*t1/(jacob*jacob);
 	//grad[2] = t1/(g11+g22+g33) - jacob*t0/((g11+g22+g33)*(g11+g22+g33));
 
@@ -6481,6 +6483,10 @@ void Octree::assign_refine_sign_quad(geoframe& geofrm, float err_tol) {
 			if (intersect_id == 1 || intersect_id == -1) {
 				eflag_on(x, y, z, level, j);
 				find_oc_id(x, y, z, level, j, intersect_id, oc_id);
+				// An edge on the volume's boundary has cells outside the grid (-1 from xyz2octcell).
+				// is_min_edge() never makes a quad for it, and get_err_grad(-1) and
+				// vtx_idx_arr_refine[-1] read and wrote out of bounds.
+				if (oc_id[1] < 0 || oc_id[2] < 0 || oc_id[3] < 0) continue;
 				
 				num_id = 0;
 				if(get_err_grad(oc_id[0]) > err_tol) num_id++;
@@ -6523,6 +6529,7 @@ void Octree::assign_refine_sign_quad(geoframe& geofrm, float err_tol) {
 				if (intersect_id == 1 || intersect_id == -1) {
 					eflag_on(x, y, z, level, j);
 					find_oc_id(x, y, z, level, j, intersect_id, oc_id);
+					if (oc_id[1] < 0 || oc_id[2] < 0 || oc_id[3] < 0) continue; // see above
 
 					num_id = 0;
 					if(vtx_idx_arr_refine[oc_id[0]] == 1) num_id++;
diff --git a/src/cvc/tests/lbie_mesher_test.cpp b/src/cvc/tests/lbie_mesher_test.cpp
index 1929122ef..25e3a48ff 100644
--- a/src/cvc/tests/lbie_mesher_test.cpp
+++ b/src/cvc/tests/lbie_mesher_test.cpp
@@ -220,6 +220,20 @@ TEST_F(LbieMesherTest, QuadDoMeshAdaptiveErrSweep) {
   expect_geoframe_valid(coarse);
 }
 
+// Octree::assign_refine_sign_quad visited every edge the surface crosses, including edges on the
+// volume's boundary, where find_oc_id gives -1 for the cells outside the grid. It then called
+// get_err_grad(-1) and read and wrote vtx_idx_arr_refine[-1]. polygonize_quad never makes a quad
+// for such an edge (is_min_edge rejects it), so those edges are now skipped. AddressSanitizer
+// flagged the read on every mesh of a surface that leaves the volume.
+TEST_F(LbieMesherTest, QuadMeshesASurfaceThatLeavesTheVolume) {
+  volume v = make_sphere_volume(ctx, 17, 1.2); // bigger than its [-1,1]^3 box
+  for (float err : {kDefaultErr, 0.01f}) {
+    LBIE::geoframe g = run_do_mesh(v, 0.0f, 0.0f, err, kDefaultErrIn, LBIE::geoframe::QUAD);
+    EXPECT_GT(g.numquads, 0) << "err=" << err;
+    expect_geoframe_valid(g);
+  }
+}
+
 TEST_F(LbieMesherTest, QuadDoMeshGeoFlow) {
   volume v = make_sphere_volume(ctx, 17, 0.5);
   LBIE::geoframe g = run_do_mesh(v, 0.0f, 0.0f, kDefaultErr, kDefaultErrIn, LBIE::geoframe::QUAD,
@@ -424,12 +438,10 @@ TEST_F(LbieMesherTest, HexaDoMeshOptimize) {
   LBIE::geoframe g = run_do_mesh(v, 0.0f, 0.0f, kDefaultErr, kDefaultErrIn, LBIE::geoframe::HEXA,
                                  LBIE::Mesher::OPTIMIZATION, 1);
   EXPECT_GT(g.numquads, 0);
-  // OPTIMIZATION is subject to the same pre-existing zero-area-division NaN
-  // bug as the geometric_flow_* improvement paths (see
-  // expect_geoframe_topology_valid): it produces finite coords on x86_64 and
-  // arm64 Debug, but a couple of NaN vertices on arm64 Release. Check only
-  // structural validity, matching the other improvement-path tests here.
-  expect_geoframe_topology_valid(g);
+  // This used to check only structure: arm64 Release left a couple of NaN vertices. Octree::getGrad
+  // read v1[3], v2[3] and v3[3] of float[3] arrays (AddressSanitizer: stack-buffer-overflow here),
+  // and v3[3] could be the uninitialized point the line search writes next. Its indices are fixed.
+  expect_geoframe_valid(g);
 }
 
 // ---------------------------------------------------------------------------