From cd71c2c1e601967029eed58516b5b626c499211a Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Thu, 27 Aug 2026 19:41:43 -0700
Subject: [PATCH 01/12] wip: shape derivatives for geometric objects

geometry_addgradient finite-differences eff_chi1inv_row with respect to an
object's centre or size, the same way the design gradient does with respect to
a design weight, and contracts against the adjoint field. No extra timestepping.

Two things differ from the design gradient. A design weight is local so it is
perturbed inside the point loop; a centre is global, so the perturbation hoists
out -- twelve passes for six parameters rather than twelve per point. And the
geometry cannot simply be mutated, because geom_box_tree caches bounding boxes
per node and per object; stale bounds would misreport which object owns a point
near the moved boundary, which is exactly where a shape derivative lives. The
tree is rebuilt around each perturbation.

The step has units here. FD_DEFAULT is a dimensionless perturbation of a weight;
a centre is a length, so the step is set as a fraction of a pixel. Subpixel
smoothing is required and refused if absent, since without it the permittivity
is a step function of position and a finite difference of that is not a
derivative.

State: correct in 2D TM, 6e-6 to 6e-4 against finite differences for centre and
size, and exactly zero for a block translated through a uniform medium. TE is
wrong by 1-38% and does not converge with resolution, so the off-diagonal
contraction is defective -- see the next commit message. Not yet fit to land.
---
 python/Makefile.am                     |   1 +
 python/adjoint/__init__.py             |   1 +
 python/adjoint/geometry_gradient.py    | 182 ++++++++++++++++++++
 python/adjoint/optimization_problem.py |  60 ++++++-
 python/geom.py                         |  39 ++++-
 python/meep.i                          |  42 +++++
 src/meepgeom.cpp                       | 226 +++++++++++++++++++++++++
 src/meepgeom.hpp                       |   6 +
 8 files changed, 552 insertions(+), 5 deletions(-)
 create mode 100644 python/adjoint/geometry_gradient.py

diff --git a/python/Makefile.am b/python/Makefile.am
index 985c3df58..2ba73a947 100644
--- a/python/Makefile.am
+++ b/python/Makefile.am
@@ -256,6 +256,7 @@ adjoint_PYTHON  = $(srcdir)/adjoint/__init__.py              \
                   $(srcdir)/adjoint/filter_source.py         \
                   $(srcdir)/adjoint/connectivity.py         \
                   $(srcdir)/adjoint/unfilter_design.py        \
+                  $(srcdir)/adjoint/geometry_gradient.py     \
                   $(srcdir)/adjoint/source_gradient.py       \
                   $(srcdir)/adjoint/wrapper.py               \
                   $(srcdir)/adjoint/utils.py
diff --git a/python/adjoint/__init__.py b/python/adjoint/__init__.py
index 16c1b3fdd..d71c614cc 100644
--- a/python/adjoint/__init__.py
+++ b/python/adjoint/__init__.py
@@ -21,6 +21,7 @@
 from .unfilter_design import *
 
 from . import source_gradient
+from . import geometry_gradient
 
 # JAX is an optional dependency; everything that needs it lives in `wrapper`.
 # Importing it also registers JAX as a way to differentiate objective functions,
diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
new file mode 100644
index 000000000..a03024a32
--- /dev/null
+++ b/python/adjoint/geometry_gradient.py
@@ -0,0 +1,182 @@
+"""Adjoint gradients with respect to a geometric object's centre and size.
+
+The design gradient finite-differences `eff_chi1inv_row` -- the subpixel
+smoothed inverse permittivity over a voxel -- with respect to a design weight.
+That function is a function of the *geometry*, so perturbing an object's centre
+or size instead gives the shape derivative from the same machinery, at no cost
+in additional timestepping. It is the same pattern as the design weights and as
+the Gaussian beam parameters: finite-difference a cheap analytic map, contract
+against the adjoint field.
+
+Two things differ from the design gradient.
+
+A design weight is local, so it can be perturbed inside the point loop. A centre
+or a size is global, so the perturbation hoists out: move the geometry once,
+visit every point, restore. Twelve passes for six parameters, not twelve
+perturbations per point.
+
+And the step has units. `utils.FD_DEFAULT` is a dimensionless perturbation of a
+weight; here it is a length, and the scale that makes sense is a fraction of a
+pixel. `DEFAULT_STEP_PIXELS` sets it relative to the grid rather than absolutely.
+"""
+
+from typing import List, Optional
+
+import numpy as np
+
+import meep as mp
+
+# Perturbation for the geometric finite difference, as a fraction of a pixel.
+# An absolute step is the wrong idea: subpixel smoothing varies over a pixel, so
+# the step has to be defined against the grid. Too small and the smoothed
+# permittivity has not changed measurably; too large and the difference samples
+# curvature rather than a derivative.
+DEFAULT_STEP_PIXELS = 0.02
+
+# Indices `geometry_addgradient` uses, matching the geom_param enumeration in
+# meepgeom.cpp.
+PARAMETER_INDEX = {
+    "center": (0, 1, 2),
+    "size": (3, 4, 5),
+}
+
+DIFFERENTIABLE_PARAMS = tuple(PARAMETER_INDEX)
+
+
+def differentiable_objects(sim: mp.Simulation) -> List:
+    """The geometric objects flagged for differentiation, with their indices."""
+    return [
+        (index, obj)
+        for index, obj in enumerate(sim.geometry)
+        if getattr(obj, "differentiable", ())
+    ]
+
+
+def object_key(obj, index: int):
+    """The key an object's gradient appears under: its name, or its position."""
+    name = getattr(obj, "name", None)
+    return name if name is not None else index
+
+
+def _check_supported(obj) -> None:
+    """Refuse the cases the diagonal-only contraction would get wrong."""
+    if not isinstance(obj, mp.Block):
+        raise NotImplementedError(
+            f"Geometry gradients are implemented for mp.Block, not "
+            f"{type(obj).__name__}. 'center' is meaningful for any object, but "
+            "'size' is not, and the contraction currently assumes axis-aligned "
+            "faces."
+        )
+    axes = (obj.e1, obj.e2, obj.e3)
+    expected = (mp.Vector3(1, 0, 0), mp.Vector3(0, 1, 0), mp.Vector3(0, 0, 1))
+    for axis, unit in zip(axes, expected):
+        if abs(axis.x - unit.x) + abs(axis.y - unit.y) + abs(axis.z - unit.z) > 1e-12:
+            raise NotImplementedError(
+                "Geometry gradients assume an axis-aligned block, so that "
+                "'size' means what the three components suggest. The "
+                "off-diagonal terms subpixel smoothing produces at a tilted "
+                "boundary *are* contracted; it is the parameterization that is "
+                "restricted, not the physics."
+            )
+
+
+def check_smoothing(sim: mp.Simulation) -> None:
+    """Refuse to differentiate a staircase.
+
+    Without subpixel smoothing the effective permittivity is a step function of
+    an object's position: shifting a boundary changes nothing until it crosses a
+    voxel edge, then changes by the full material contrast. A finite difference
+    of that is zero almost everywhere and a spike in a few voxels, and the sum
+    is not a derivative. It would look like a plausible number rather than an
+    error, so refuse instead.
+    """
+    if not getattr(sim, "eps_averaging", True):
+        raise ValueError(
+            "Geometry gradients need subpixel smoothing, but this simulation "
+            "has eps_averaging=False. Without it the permittivity is a step "
+            "function of an object's position and the derivative does not "
+            "exist at the grid scale."
+        )
+
+
+def install_geometry_monitors(
+    sim: mp.Simulation,
+    objects,
+    frequencies,
+    decimation_factor: int = 0,
+) -> List[List[mp.DftFields]]:
+    """DFT monitors over each differentiable object's support.
+
+    Unlike a design region, an ordinary object has no monitor of its own, so the
+    fields it needs have to be recorded in both the forward and adjoint runs.
+    """
+    from . import utils
+
+    monitors = []
+    for _, obj in objects:
+        volume = sim._fit_volume_to_simulation(
+            mp.Volume(center=obj.center, size=obj.size)
+        )
+        monitors.append(
+            [
+                sim.add_dft_fields(
+                    [component],
+                    frequencies,
+                    where=volume,
+                    yee_grid=True,
+                    decimation_factor=decimation_factor,
+                    persist=True,
+                )
+                for component in utils._compute_components(sim)
+            ]
+        )
+    return monitors
+
+
+def gradient(
+    sim: mp.Simulation,
+    obj,
+    object_index: int,
+    forward_fields: List[mp.DftFields],
+    adjoint_fields: List[mp.DftFields],
+    frequencies,
+    step: Optional[float] = None,
+) -> dict:
+    """dJ/d(parameter) for one object, keyed by the names it declared."""
+    _check_supported(obj)
+    check_smoothing(sim)
+
+    names = [n for n in obj.differentiable]
+    indices = []
+    for name in names:
+        indices.extend(PARAMETER_INDEX[name])
+    indices = np.asarray(indices, dtype=np.intc)
+
+    frequencies = np.asarray(frequencies, dtype=np.float64)
+    out = np.zeros((frequencies.size, indices.size), dtype=np.float64)
+
+    if step is None:
+        step = DEFAULT_STEP_PIXELS / sim.resolution
+
+    mp._get_geometry_gradient(
+        out,
+        1.0,
+        adjoint_fields[0].swigobj,
+        adjoint_fields[1].swigobj,
+        adjoint_fields[2].swigobj,
+        forward_fields[0].swigobj,
+        forward_fields[1].swigobj,
+        forward_fields[2].swigobj,
+        sim.gv,
+        frequencies,
+        sim.geps,
+        object_index,
+        indices,
+        step,
+    )
+
+    grads = {}
+    for position, name in enumerate(names):
+        columns = slice(3 * position, 3 * position + 3)
+        grads[name] = np.squeeze(out[:, columns])
+    return grads
diff --git a/python/adjoint/optimization_problem.py b/python/adjoint/optimization_problem.py
index dc6a38097..36b387e25 100644
--- a/python/adjoint/optimization_problem.py
+++ b/python/adjoint/optimization_problem.py
@@ -6,6 +6,7 @@
 
 from . import LDOS, DesignRegion, utils, ObjectiveQuantity
 from . import source_gradient
+from . import geometry_gradient
 
 
 class OptimizationProblem:
@@ -154,6 +155,12 @@ def __init__(
         self.differentiable_sources = source_gradient.differentiable_sources(self.sim)
         self.source_gradient = {}
 
+        # Geometric objects flagged with `differentiable=[...]`: the shape
+        # derivative of the structure itself, as opposed to the density inside
+        # a design region.
+        self.differentiable_objects = geometry_gradient.differentiable_objects(self.sim)
+        self.geometry_gradient = {}
+
         # The optimizer has three allowable states : "INIT", "FWD", and "ADJ".
         #    INIT - The optimizer is initialized and ready to run a forward simulation
         #    FWD  - The optimizer has already run a forward simulation
@@ -227,10 +234,14 @@ def __call__(
                     f"Incorrect solver state detected: {self.current_state}"
                 )
 
-        if self.differentiable_sources:
-            # Only change the return shape when the user asked for source
-            # gradients; without a flagged source this is exactly as before.
-            return self.f0, {"design": self.gradient, **self.source_gradient}
+        if self.differentiable_sources or self.differentiable_objects:
+            # Only change the return shape when something extra was flagged;
+            # otherwise this is exactly as before.
+            return self.f0, {
+                "design": self.gradient,
+                **self.source_gradient,
+                **self.geometry_gradient,
+            }
 
         return self.f0, self.gradient
 
@@ -270,6 +281,14 @@ def prepare_forward_run(self):
         self.forward_design_region_monitors = utils.install_design_region_monitors(
             self.sim, self.design_regions, self.frequencies, self.decimation_factor
         )
+        # an ordinary object has no monitor of its own, so its fields have to be
+        # recorded in both runs
+        self.forward_geometry_monitors = geometry_gradient.install_geometry_monitors(
+            self.sim,
+            self.differentiable_objects,
+            self.frequencies,
+            self.decimation_factor,
+        )
 
     def forward_run(self):
         # set up monitors
@@ -361,6 +380,7 @@ def adjoint_run(self):
 
         self.adjoint_design_region_monitors = []
         self.adjoint_source_monitors = []
+        self.adjoint_geometry_monitors = []
         for ar in range(len(self.objective_functions)):
             # Reset the fields
             self.sim.restart_fields()
@@ -381,6 +401,15 @@ def adjoint_run(self):
 
             # register a monitor over each differentiable source's support; the
             # adjoint field there is the gradient with respect to its currents
+            self.adjoint_geometry_monitors.append(
+                geometry_gradient.install_geometry_monitors(
+                    self.sim,
+                    self.differentiable_objects,
+                    self.frequencies,
+                    self.decimation_factor,
+                )
+            )
+
             self.adjoint_source_monitors.append(
                 source_gradient.install_source_gradient_monitors(
                     self.sim,
@@ -498,8 +527,31 @@ def _source_gradient_for(self, src, monitors, scale):
                 )
         return grads
 
+    def calculate_geometry_gradient(self):
+        """Shape derivatives for each flagged geometric object."""
+        if not self.differentiable_objects:
+            return {}
+        out = {}
+        for ar in range(len(self.objective_functions)):
+            for oi, (object_index, obj) in enumerate(self.differentiable_objects):
+                grads = geometry_gradient.gradient(
+                    self.sim,
+                    obj,
+                    object_index,
+                    self.forward_geometry_monitors[oi],
+                    self.adjoint_geometry_monitors[ar][oi],
+                    self.frequencies,
+                )
+                key = geometry_gradient.object_key(obj, object_index)
+                if len(self.objective_functions) == 1:
+                    out[key] = grads
+                else:
+                    out.setdefault(key, []).append(grads)
+        return out
+
     def calculate_gradient(self):
         self.source_gradient = self.calculate_source_gradient()
+        self.geometry_gradient = self.calculate_geometry_gradient()
 
         # Iterate through all design regions and calculate gradient
         self.gradient = [
diff --git a/python/geom.py b/python/geom.py
index b1207c75b..b15da1ce2 100755
--- a/python/geom.py
+++ b/python/geom.py
@@ -1033,6 +1033,35 @@ def __init__(
         self.pumping_rate = pumping_rate
 
 
+# Parameters a geometric object can be differentiated with respect to.
+# Unlike a source, where `center` and `size` move which grid points are driven
+# and so fall outside the formulation, epsilon is a function of position and
+# moving an object is exactly what a derivative with respect to position means.
+_DIFFERENTIABLE_GEOMETRY = ("center", "size")
+
+
+def _validate_differentiable_geometry(differentiable):
+    """Check and normalize a geometric object's `differentiable` argument."""
+    if differentiable is None:
+        return ()
+    if isinstance(differentiable, str):
+        raise ValueError(
+            "`differentiable` takes a list of parameter names, not a bare "
+            f"string; use ['{differentiable}'] instead."
+        )
+    names = list(differentiable)
+    for name in names:
+        if name not in _DIFFERENTIABLE_GEOMETRY:
+            raise ValueError(
+                f"'{name}' is not a differentiable parameter of a geometric "
+                f"object. Valid choices are: "
+                f"{', '.join(_DIFFERENTIABLE_GEOMETRY)}."
+            )
+    if len(set(names)) != len(names):
+        raise ValueError(f"`differentiable` contains duplicate names: {names}")
+    return tuple(names)
+
+
 class GeometricObject:
     """
     This class, and its descendants, are used to specify the solid geometric objects that
@@ -1088,7 +1117,13 @@ class GeometricObject:
     """
 
     def __init__(
-        self, material=Medium(), center=Vector3(), epsilon_func=None, label=None
+        self,
+        material=Medium(),
+        center=Vector3(),
+        epsilon_func=None,
+        label=None,
+        differentiable=None,
+        name=None,
     ):
         """
         Construct a `GeometricObject`.
@@ -1119,6 +1154,8 @@ def __init__(
 
         self.label = label
         self.material = material
+        self.name = name
+        self.differentiable = _validate_differentiable_geometry(differentiable)
         self.center = Vector3(*center)
 
     def __contains__(self, point):
diff --git a/python/meep.i b/python/meep.i
index 51fa5aa3d..8b331eb97 100644
--- a/python/meep.i
+++ b/python/meep.i
@@ -652,6 +652,7 @@ void _get_eigenmode(meep::fields *f, double frequency, meep::direction d, const
 %feature("nothreadallow") py_do_harminv;
 %feature("nothreadallow") _get_array_slice_dimensions;
 %feature("nothreadallow") _get_gradient;
+%feature("nothreadallow") _get_geometry_gradient;
 %feature("nothreadallow") _get_dft_array;
 
 %numpy_typemaps(std::complex<double>, NPY_CDOUBLE, int);
@@ -866,6 +867,47 @@ meep::volume_list *make_volume_list(const meep::volume &v, int c,
 // typemaps needed for material grid
 //--------------------------------------------------
 
+%inline %{
+void _get_geometry_gradient(PyObject *grad, double scalegrad,
+                   meep::dft_fields *fields_a_0, meep::dft_fields *fields_a_1, meep::dft_fields *fields_a_2,
+                   meep::dft_fields *fields_f_0, meep::dft_fields *fields_f_1, meep::dft_fields *fields_f_2,
+                   meep::grid_volume *grid_volume, PyObject *frequencies,
+                   meep_geom::geom_epsilon *geps, int object_index, PyObject *params,
+                   double fd_step) {
+
+    PyArrayObject *pao_grad = (PyArrayObject *)grad;
+    if (!PyArray_Check(pao_grad)) meep::abort("grad parameter must be numpy array.");
+    if (!PyArray_ISCARRAY(pao_grad)) meep::abort("Numpy grad array must be C-style contiguous.");
+    if (PyArray_NDIM(pao_grad) != 2) meep::abort("Numpy grad array must have 2 dimensions.");
+    double *grad_c = (double *)PyArray_DATA(pao_grad);
+    npy_intp nparams = PyArray_DIMS(pao_grad)[1];
+
+    PyArrayObject *pao_params = (PyArrayObject *)params;
+    if (!PyArray_Check(pao_params)) meep::abort("params must be a numpy array.");
+    if (!PyArray_ISCARRAY(pao_params)) meep::abort("Numpy params array must be C-style contiguous.");
+    if (PyArray_DIMS(pao_params)[0] != nparams)
+        meep::abort("params has %td entries but grad is allocated for %td.",
+                    PyArray_DIMS(pao_params)[0], nparams);
+    int *params_c = (int *)PyArray_DATA(pao_params);
+
+    std::vector<meep::dft_fields *> adjoint_fields = {fields_a_0,fields_a_1,fields_a_2};
+    std::vector<meep::dft_fields *> forward_fields = {fields_f_0,fields_f_1,fields_f_2};
+
+    PyArrayObject *pao_freqs = (PyArrayObject *)frequencies;
+    if (!PyArray_Check(pao_freqs)) meep::abort("frequencies parameter must be numpy array.");
+    if (!PyArray_ISCARRAY(pao_freqs)) meep::abort("Numpy frequencies array must be C-style contiguous.");
+    double *frequencies_c = (double *)PyArray_DATA(pao_freqs);
+    npy_intp nf = PyArray_DIMS(pao_freqs)[0];
+    if (PyArray_DIMS(pao_grad)[0] != nf)
+        meep::abort("Numpy grad array is allocated for %td frequencies; it should be allocated for %td.",
+                    PyArray_DIMS(pao_grad)[0], nf);
+
+    meep_geom::geometry_addgradient(grad_c, nparams, nf, adjoint_fields, forward_fields,
+                                    frequencies_c, scalegrad, *grid_volume, geps,
+                                    object_index, params_c, fd_step);
+}
+%}
+
 %inline %{
 void _get_gradient(PyObject *grad, double scalegrad,
                    meep::dft_fields *fields_a_0, meep::dft_fields *fields_a_1, meep::dft_fields *fields_a_2,
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index dce99195e..23973a191 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -2940,6 +2940,232 @@ static std::complex<meep::realnum> forward_dft_value(const meep::dft_chunk *ch,
     if (i >= 0 && (size_t)i < ch->N) return ch->dft[nf * i + f_i];
   }
   return 0;
+
+/* ------------------------------------------------------------------ */
+/* Gradients with respect to a geometric object's centre and size.      */
+/* ------------------------------------------------------------------ */
+
+/* The design gradient finite-differences eff_chi1inv_row with respect to a
+   design weight. eff_chi1inv_row is a function of the *geometry*, so
+   perturbing an object's centre or size instead gives dA/d(parameter) from the
+   same machinery, with no additional timestepping.
+
+   One structural difference. A design weight is local, so the design gradient
+   perturbs it inside the point loop. A centre or a size is global, so the
+   perturbation hoists out: the geometry is moved once, every point is visited,
+   and the geometry is restored. That is twelve passes for six parameters, not
+   twelve perturbations per point.
+
+   The perturbation cannot simply mutate the object, because geom_box_tree
+   caches bounding boxes (`geom_box b, b1, b2` per node, plus one per object).
+   Stale bounds would misreport which object owns a point near the moved
+   boundary -- exactly where a shape derivative lives -- so the tree is rebuilt
+   around each perturbation. */
+
+// Which scalar of an object is being perturbed.
+enum geom_param {
+  GEOM_CENTER_X = 0,
+  GEOM_CENTER_Y,
+  GEOM_CENTER_Z,
+  GEOM_SIZE_X,
+  GEOM_SIZE_Y,
+  GEOM_SIZE_Z
+};
+
+static double *geom_param_slot(geometric_object *o, int which) {
+  switch (which) {
+    case GEOM_CENTER_X: return &o->center.x;
+    case GEOM_CENTER_Y: return &o->center.y;
+    case GEOM_CENTER_Z: return &o->center.z;
+    default: break;
+  }
+  if (o->which_subclass != geometric_object::BLOCK) return NULL;
+  switch (which) {
+    case GEOM_SIZE_X: return &o->subclass.block_data->size.x;
+    case GEOM_SIZE_Y: return &o->subclass.block_data->size.y;
+    case GEOM_SIZE_Z: return &o->subclass.block_data->size.z;
+    default: return NULL;
+  }
+}
+
+/* Rebuild the box tree after the geometry underneath it has moved. */
+static void geom_rebuild_tree(geom_epsilon *geps, const meep::grid_volume &gv) {
+  if (geps->restricted_tree && geps->restricted_tree != geps->geometry_tree)
+    destroy_geom_box_tree(geps->restricted_tree);
+  destroy_geom_box_tree(geps->geometry_tree);
+  geom_fix_object_list(geps->geometry);
+  geom_box box = gv2box(gv.surroundings());
+  geps->geometry_tree = create_geom_box_tree0(geps->geometry, box);
+  geps->restricted_tree = geps->geometry_tree;
+}
+
+/* One point's contribution: the smoothed chi1inv row at `r`, contracted with
+   the adjoint and forward fields. */
+static std::complex<double> geom_chi1inv_term(geom_epsilon *geps, const meep::grid_volume &gv,
+                                              meep::component adjoint_c, int dir_idx,
+                                              const meep::vec &r, double scale,
+                                              std::complex<meep::realnum> adj,
+                                              std::complex<meep::realnum> fwd) {
+  meep::volume voxel(r);
+  LOOP_OVER_DIRECTIONS(gv.dim, d) {
+    voxel.set_direction_min(d, r.in_direction(d) - 0.5 * gv.inva);
+    voxel.set_direction_max(d, r.in_direction(d) + 0.5 * gv.inva);
+  }
+  double row[3];
+  geps->eff_chi1inv_row(adjoint_c, row, voxel, geps->tol, geps->maxeval);
+  return scale * row[dir_idx] * std::complex<double>(double(adj.real()), double(adj.imag())) *
+         std::complex<double>(double(fwd.real()), double(fwd.imag()));
+}
+
+/* Contract the adjoint and forward fields against the *current* geometry,
+   accumulating one number per frequency. Called twice per parameter, on either
+   side of the perturbation. */
+static void geom_contract(std::complex<double> *out, size_t nf,
+                          std::vector<meep::dft_chunk *> *adjoint_chunks,
+                          std::vector<meep::dft_chunk *> *forward_chunks, double scalegrad,
+                          meep::grid_volume &gv, geom_epsilon *geps) {
+  for (size_t f_i = 0; f_i < nf; f_i++) {
+    for (int ci_adjoint = 0; ci_adjoint < 3; ci_adjoint++) {
+      int num_chunks = adjoint_chunks[ci_adjoint].size();
+      if (num_chunks == 0) continue;
+
+      for (int cur = 0; cur < num_chunks; cur++) {
+        meep::dft_chunk *adj_chunk = adjoint_chunks[ci_adjoint][cur];
+        meep::component adjoint_c = adj_chunk->c;
+        meep::grid_volume gv_adj = gv.subvolume(adj_chunk->is, adj_chunk->ie, adjoint_c);
+
+        for (int ci_forward = 0; ci_forward < 3; ci_forward++) {
+          size_t num_f = forward_chunks[ci_forward].size();
+          if ((num_f == 0) || ((size_t)cur >= num_f)) continue;
+          meep::dft_chunk *fwd_chunk = forward_chunks[ci_forward][cur];
+          meep::component forward_c = fwd_chunk->c;
+          meep::grid_volume gv_fwd = gv.subvolume(fwd_chunk->is, fwd_chunk->ie, forward_c);
+
+          int dir_idx;
+          switch (meep::component_direction(forward_c)) {
+            case meep::X:
+            case meep::R: dir_idx = 0; break;
+            case meep::Y:
+            case meep::P: dir_idx = 1; break;
+            case meep::Z: dir_idx = 2; break;
+            default: continue;
+          }
+
+          LOOP_OVER_IVECS(gv_adj, adj_chunk->is_old, adj_chunk->ie_old, idx_adj) {
+            IVEC_LOOP_ILOC(gv_adj, ip);
+            IVEC_LOOP_LOC(gv_adj, p);
+            std::complex<meep::realnum> adj = adj_chunk->dft[nf * idx_adj + f_i];
+            if (adj == 0.0) continue;
+
+            if (forward_c == adjoint_c) {
+              std::complex<meep::realnum> fwd = fwd_chunk->dft[nf * idx_adj + f_i];
+              double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p.r() : 1;
+              out[f_i] += geom_chi1inv_term(geps, gv, adjoint_c, dir_idx, p, scalegrad * cyl_scale,
+                                            adj, fwd);
+            }
+            else {
+              /* Subpixel smoothing is second order because it builds an
+                 effective *tensor*: harmonic averaging for the field component
+                 normal to the interface, arithmetic for the tangential ones.
+                 The off-diagonal entries that produces are non-zero wherever
+                 the interface normal is not along an axis -- which for a shape
+                 derivative is the whole of the signal, since the derivative
+                 lives entirely at the boundary.
+
+                 Contracting them needs the forward field of a different
+                 component, which sits at a different Yee location, so it is
+                 restricted to the two epsilon nodes between the pair and
+                 interpolated there. This mirrors material_grids_addgradient. */
+              meep::ivec fwd_p = ip + gv.iyee_shift(forward_c) - gv.iyee_shift(adjoint_c);
+              meep::ivec unit_a = unit_ivec(gv.dim, component_direction(adjoint_c));
+              meep::ivec unit_f = unit_ivec(gv.dim, component_direction(forward_c));
+              meep::ivec fwd_pa = fwd_p + unit_a * 2;
+              meep::ivec fwd_pf = fwd_p - unit_f * 2;
+              meep::ivec fwd_paf = fwd_p + unit_a * 2 - unit_f * 2;
+
+              meep::ivec fwd_pl[2] = {fwd_p, fwd_pa};
+              meep::ivec fwd_pr[2] = {fwd_pf, fwd_paf};
+              meep::ivec ieps[2] = {(fwd_p + fwd_pf) / 2, (fwd_pa + fwd_paf) / 2};
+
+              for (int node = 0; node < 2; node++) {
+                ptrdiff_t i1 = gv_fwd.index(forward_c, fwd_pl[node]);
+                ptrdiff_t i2 = gv_fwd.index(forward_c, fwd_pr[node]);
+                std::complex<meep::realnum> fwd1 =
+                    ((i1 >= fwd_chunk->N) || (i1 < 0)) ? 0 : fwd_chunk->dft[nf * i1 + f_i];
+                std::complex<meep::realnum> fwd2 =
+                    ((i2 >= fwd_chunk->N) || (i2 < 0)) ? 0 : fwd_chunk->dft[nf * i2 + f_i];
+                std::complex<meep::realnum> fwd_avg =
+                    std::complex<meep::realnum>(0.5, 0) * (fwd1 + fwd2);
+                meep::vec eps1 = gv[ieps[node]];
+                double cyl_scale = (gv.dim == meep::Dcyl) ? eps1.r() : 1;
+                out[f_i] +=
+                    geom_chi1inv_term(geps, gv, adjoint_c, dir_idx, eps1, scalegrad * cyl_scale,
+                                      std::complex<meep::realnum>(0.5, 0) * adj, fwd_avg);
+              }
+            }
+          }
+        }
+      }
+    }
+  }
+}
+
+/* dJ/d(parameter) for one geometric object.
+
+   `v` receives nf * nparams doubles, frequency-major. `params` names which
+   scalars to differentiate, using the geom_param enumeration.
+
+   Both the diagonal and the off-diagonal terms of the smoothed tensor are
+   contracted; see geom_contract for why the second matters here more than it
+   does for a density gradient. */
+void geometry_addgradient(double *v, size_t nparams, size_t nf,
+                          std::vector<meep::dft_fields *> fields_a,
+                          std::vector<meep::dft_fields *> fields_f, double *frequencies,
+                          double scalegrad, meep::grid_volume &gv, geom_epsilon *geps,
+                          int object_index, int *params, double du) {
+  (void)frequencies;
+  if (object_index < 0 || object_index >= geps->geometry.num_items)
+    meep::abort("geometry_addgradient: object index %d out of range (%d objects)", object_index,
+                geps->geometry.num_items);
+  geometric_object *obj = &geps->geometry.items[object_index];
+
+  std::vector<meep::dft_chunk *> adjoint_chunks[3], forward_chunks[3];
+  for (int i = 0; i < 3; i++) {
+    for (meep::dft_chunk *c = fields_a[i]->chunks; c; c = c->next_in_dft)
+      adjoint_chunks[i].push_back(c);
+    for (meep::dft_chunk *c = fields_f[i]->chunks; c; c = c->next_in_dft)
+      forward_chunks[i].push_back(c);
+  }
+
+  std::vector<std::complex<double> > plus(nf), minus(nf);
+  std::vector<double> local(nf * nparams, 0.0);
+
+  for (size_t ip = 0; ip < nparams; ip++) {
+    double *slot = geom_param_slot(obj, params[ip]);
+    if (!slot)
+      meep::abort("geometry_addgradient: parameter %d is not available on this object", params[ip]);
+    double original = *slot;
+
+    std::fill(minus.begin(), minus.end(), std::complex<double>(0, 0));
+    std::fill(plus.begin(), plus.end(), std::complex<double>(0, 0));
+
+    *slot = original - du;
+    geom_rebuild_tree(geps, gv);
+    geom_contract(minus.data(), nf, adjoint_chunks, forward_chunks, scalegrad, gv, geps);
+
+    *slot = original + du;
+    geom_rebuild_tree(geps, gv);
+    geom_contract(plus.data(), nf, adjoint_chunks, forward_chunks, scalegrad, gv, geps);
+
+    *slot = original;
+    geom_rebuild_tree(geps, gv);
+
+    // (row_1 - row_2) / 2du, matching get_material_gradient's sign convention
+    for (size_t f_i = 0; f_i < nf; f_i++)
+      local[nparams * f_i + ip] = std::real(minus[f_i] - plus[f_i]) / (2 * du);
+  }
+
+  meep::sum_to_all(local.data(), v, int(nf * nparams));
 }
 
 void material_grids_addgradient(double *v, size_t ng, size_t nf,
diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp
index 004954209..ce60e1577 100644
--- a/src/meepgeom.hpp
+++ b/src/meepgeom.hpp
@@ -294,6 +294,12 @@ meep::vec material_grid_grad(vector3 p, material_data *md, const geometric_objec
 double matgrid_val(vector3 p, geom_box_tree tp, int oi, material_data *md);
 double material_grid_val(vector3 p, material_data *md);
 geom_box_tree calculate_tree(const meep::volume &v, geometric_object_list g);
+void geometry_addgradient(double *v, size_t nparams, size_t nf,
+                          std::vector<meep::dft_fields *> fields_a,
+                          std::vector<meep::dft_fields *> fields_f, double *frequencies,
+                          double scalegrad, meep::grid_volume &gv, geom_epsilon *geps,
+                          int object_index, int *params, double du);
+
 void material_grids_addgradient(double *v, size_t ng, size_t nf,
                                 std::vector<meep::dft_fields *> fields_a,
                                 std::vector<meep::dft_fields *> fields_f, double *frequencies,

From 1fc4fc06d8f74892b5a64dc0753677804df19eb2 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Thu, 27 Aug 2026 20:25:17 -0700
Subject: [PATCH 02/12] wip: narrow the TE failure in the shape derivative

Three hypotheses eliminated, none of them the cause.

The off-diagonal tensor terms are dead code: disabling them is bit-identical in
both polarizations. Meep's own design gradient is accurate in TE (4.8e-5)
without ever taking that branch either, so second-order smoothing's
off-diagonal entries are not what TE needs here. The contraction is kept
because it is correct, not because it is load-bearing.

The design gradient itself is not at fault: same geometry, same polarization,
differentiating the density instead of the boundary gives 4.8e-5.

Monitor padding is not the cause. The support of d(epsilon)/d(parameter) is the
shell of voxels a boundary sweeps through, which extends outside the object, so
padding by two pixels is right in principle -- it grows the monitor from 437 to
621 points -- but the added points contribute exactly zero and TE is unchanged.
Kept as correct-in-principle, with no evidence it is needed.

TM remains accurate (6e-6 to 6e-4, and exactly zero for a block translated
through a uniform medium). TE is 1-38% and diverges with resolution.
---
 python/adjoint/geometry_gradient.py | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
index a03024a32..0aa6e46b0 100644
--- a/python/adjoint/geometry_gradient.py
+++ b/python/adjoint/geometry_gradient.py
@@ -33,6 +33,11 @@
 # curvature rather than a derivative.
 DEFAULT_STEP_PIXELS = 0.02
 
+# How far past an object the gradient monitor reaches, in pixels. The support of
+# d(epsilon)/d(parameter) is the shell of voxels the boundary sweeps through, so
+# it extends outside the object itself.
+MONITOR_PAD_PIXELS = 2.0
+
 # Indices `geometry_addgradient` uses, matching the geom_param enumeration in
 # meepgeom.cpp.
 PARAMETER_INDEX = {
@@ -114,8 +119,20 @@ def install_geometry_monitors(
 
     monitors = []
     for _, obj in objects:
+        # Pad beyond the object. Moving a boundary changes the smoothed
+        # permittivity in every voxel the boundary passes through, and that
+        # includes voxels whose centres lie *outside* the original extent. A
+        # monitor covering only the object misses them, which costs little
+        # where the field is continuous across the boundary and a great deal
+        # where it is not -- the normal E field jumps by the index contrast.
+        pad = MONITOR_PAD_PIXELS / sim.resolution
+        padded = mp.Vector3(
+            obj.size.x + 2 * pad if obj.size.x else 0.0,
+            obj.size.y + 2 * pad if obj.size.y else 0.0,
+            obj.size.z + 2 * pad if obj.size.z else 0.0,
+        )
         volume = sim._fit_volume_to_simulation(
-            mp.Volume(center=obj.center, size=obj.size)
+            mp.Volume(center=obj.center, size=padded)
         )
         monitors.append(
             [

From 7d8bc385acfd6c26d426474b9e3ee76f70f522b7 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 08:16:11 -0700
Subject: [PATCH 03/12] wip: the shape derivative's real failure mode is a kink
 at voxel edges

The step sweep from the design doc, which I had not run, explains the 2D
failures. There is no plateau when a block face sits on a voxel edge: the
reported gradient sweeps monotonically from -0.19 to +0.12 as the step goes
from 0.2 px to 0.001 px, passing through the true value without settling.
Moving the same face half a pixel makes the sweep converge -- 10.1, 9.6, 6.5,
3.4, 1.3, 0.72 percent over that same range.

A voxel's filling fraction is piecewise linear in the position of the boundary
crossing it, kinked wherever the boundary reaches a voxel edge, so epsilon and
the objective are C0 but not C1 in an object's position. On an edge the left
and right derivatives differ and no derivative exists; a central difference
returns a step-dependent mixture of the two. This is easy to hit by accident,
since round sizes at round resolutions land on edges -- a 0.8 block at
resolution 20 has faces exactly 8 pixels from its centre. Now warned about.

It also means the one case that looked right was a false positive. 2D TM had
its faces on edges too, but with a single field component and a symmetric
geometry the two neighbouring voxels carried near-identical fields and the
mixture cancelled.

The default step was also far too large: 0.02 px costs 3.4% on its own. Now
0.002 px.

3D is not explained by either. A smaller step moves centre.x from 38% to 17%
but it plateaus near 15%, so there is a second cause there.
---
 python/adjoint/geometry_gradient.py | 55 ++++++++++++++++++++++++++---
 1 file changed, 51 insertions(+), 4 deletions(-)

diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
index 0aa6e46b0..6f019b9c3 100644
--- a/python/adjoint/geometry_gradient.py
+++ b/python/adjoint/geometry_gradient.py
@@ -20,6 +20,7 @@
 pixel. `DEFAULT_STEP_PIXELS` sets it relative to the grid rather than absolutely.
 """
 
+import warnings
 from typing import List, Optional
 
 import numpy as np
@@ -28,10 +29,13 @@
 
 # Perturbation for the geometric finite difference, as a fraction of a pixel.
 # An absolute step is the wrong idea: subpixel smoothing varies over a pixel, so
-# the step has to be defined against the grid. Too small and the smoothed
-# permittivity has not changed measurably; too large and the difference samples
-# curvature rather than a derivative.
-DEFAULT_STEP_PIXELS = 0.02
+# the step has to be defined against the grid.
+#
+# Measured on a block whose faces sit mid-voxel, against a finite difference of
+# the objective: 0.2 px is 10% out, 0.05 px is 6.5%, 0.02 px is 3.4%, 0.005 px
+# is 1.3% and 0.001 px is 0.7%, still improving. The error is truncation, not
+# roundoff, over that whole range, so the step wants to be small.
+DEFAULT_STEP_PIXELS = 0.002
 
 # How far past an object the gradient monitor reaches, in pixels. The support of
 # d(epsilon)/d(parameter) is the shell of voxels the boundary sweeps through, so
@@ -85,6 +89,48 @@ def _check_supported(obj) -> None:
             )
 
 
+def check_faces_off_voxel_edges(sim: mp.Simulation, obj) -> None:
+    """Warn when a face sits exactly on a voxel edge, where no derivative exists.
+
+    A voxel's filling fraction is piecewise linear in the position of the
+    boundary crossing it, with a kink each time the boundary reaches a voxel
+    edge. The smoothed permittivity -- and so the objective -- is therefore C0
+    but not C1 in an object's position, and *on* an edge the left and right
+    derivatives differ.
+
+    A central difference straddles that kink and returns a weighted mixture of
+    the two, with weights that depend on the step. Measured on such a face, the
+    reported gradient swept monotonically from -0.19 to +0.12 as the step went
+    from 0.2 px to 0.001 px, passing through the true value without settling.
+    Moving the same face half a pixel made the sweep converge to 0.7%.
+
+    This is easy to hit by accident, because round sizes at round resolutions
+    land on edges: a 0.8-wide block at resolution 20 has faces exactly 8 pixels
+    from its centre.
+    """
+    dx = 1.0 / sim.resolution
+    size = (obj.size.x, obj.size.y, obj.size.z)
+    center = (obj.center.x, obj.center.y, obj.center.z)
+    for axis, letter in enumerate("xyz"):
+        if size[axis] == 0:
+            continue
+        for face in (center[axis] - size[axis] / 2, center[axis] + size[axis] / 2):
+            offset = abs((face / dx) - round(face / dx))
+            if offset < 1e-6:
+                warnings.warn(
+                    f"The {letter} face of "
+                    f"{getattr(obj, 'name', None) or 'this object'} at "
+                    f"{face:.6g} lies on a voxel edge, where the smoothed "
+                    "permittivity has a kink and the derivative with respect "
+                    "to position is one-sided. The reported gradient will "
+                    "depend on the finite-difference step rather than "
+                    f"converging. Offset the geometry by about {dx / 2:.6g} "
+                    "along " + letter + ", or change the resolution.",
+                    RuntimeWarning,
+                    stacklevel=3,
+                )
+
+
 def check_smoothing(sim: mp.Simulation) -> None:
     """Refuse to differentiate a staircase.
 
@@ -162,6 +208,7 @@ def gradient(
     """dJ/d(parameter) for one object, keyed by the names it declared."""
     _check_supported(obj)
     check_smoothing(sim)
+    check_faces_off_voxel_edges(sim, obj)
 
     names = [n for n in obj.differentiable]
     indices = []

From 1b1ae8769761ecb83164cf2eff0b83ef687ee8d0 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 09:52:41 -0700
Subject: [PATCH 04/12] wip: a shape derivative is one order worse than the
 fields

Measured rather than assumed. With smoothing off, the objective is piecewise
constant in an object's position, swinging 77% across one pixel. With it on the
objective is continuous and the swing falls to 14%, so smoothing does deliver
the smooth position dependence it promises.

The residual oscillates with a period of exactly one pixel, which identifies it
as a function of the sub-pixel phase rather than anything physical. Writing
J(p) = J_exact(p) + Delta^2 E(p/Delta) and differentiating turns the Delta^2
into Delta, because d/dp of E(p/Delta) carries a 1/Delta. So second-order
smoothing gives a first-order-accurate shape derivative.

Confirmed: the per-pixel swing is 9.9, 4.0 and 1.9 percent at resolutions 20,
40 and 80 -- halving per doubling, not quartering.

This corrects two earlier claims of mine. The kink story does not survive
arithmetic: a central difference across a kink returns the average of the
one-sided slopes and is step-independent, which is not what was measured. And
the 2D TM agreement at 6e-6 was not verification of anything -- it was a
symmetric configuration in which the artifact cancelled.

The implementation should therefore be tested against a finite difference of
the discrete objective, which is what it actually computes. Agreement with the
physically intended derivative is separately limited to O(Delta).
---
 python/adjoint/geometry_gradient.py | 32 +++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)

diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
index 6f019b9c3..0176676b8 100644
--- a/python/adjoint/geometry_gradient.py
+++ b/python/adjoint/geometry_gradient.py
@@ -18,6 +18,38 @@
 And the step has units. `utils.FD_DEFAULT` is a dimensionless perturbation of a
 weight; here it is a length, and the scale that makes sense is a fraction of a
 pixel. `DEFAULT_STEP_PIXELS` sets it relative to the grid rather than absolutely.
+
+Accuracy: a shape derivative is one order worse than the fields
+------------------------------------------------------------------
+Subpixel smoothing makes the permittivity -- and so the objective -- a
+second-order accurate, *continuous* function of an object's position. Measured
+on a block of index 2.5 at resolution 20, turning smoothing off makes the
+objective piecewise constant in position, swinging 77% across a single pixel;
+turning it on makes it continuous, swinging 14%.
+
+But second-order accuracy in the value is only first-order accuracy in the
+derivative. Writing the discretization error as a function of the sub-pixel
+phase,
+
+    J(p) = J_exact(p) + Delta^2 E(p / Delta)
+
+with E an O(1) oscillatory function -- which is why the residual has a period of
+exactly one pixel -- differentiating gives
+
+    dJ/dp = dJ_exact/dp + Delta^2 (1/Delta) E'(p / Delta)
+          = dJ_exact/dp + O(Delta)
+
+Measured: the per-pixel swing falls 9.9% -> 4.0% -> 1.9% at resolutions
+20 -> 40 -> 80, halving with each doubling rather than quartering. This is
+inherent to differentiating a smoothed staircase and cannot be recovered by
+improving the adjoint, because the error is in the function being
+differentiated, not in how it is differentiated.
+
+Two consequences for testing. The adjoint should be checked against a finite
+difference of the *discrete* objective, since both differentiate the same
+smoothed function and should agree closely; that is a statement about the
+implementation. Agreement with the physically intended derivative is limited to
+O(Delta) and is a statement about the discretization.
 """
 
 import warnings

From 5a8e2780923381dd463e1b379c65971f1e6a1868 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 10:32:47 -0700
Subject: [PATCH 05/12] wip: analytic dA/dp for geometry, replacing the
 geometry finite difference

The old approach perturbed the geometry and re-differenced eff_chi1inv_row.
That cannot work: the fill fraction comes from box_overlap_with_object, an
adaptive quadrature, so differencing it amplifies its tolerance by 1/step.
Measured, there was no stable step regime at all -- the reported gradient swept
monotonically across two decades of step without settling.

The smoothed tensor depends on geometry only through the fill fraction and the
interface normal, so

    d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp

d(fill)/dp is analytic and separable, because a block is an intersection of
three slabs and the pixel overlap factorizes into clamped one-dimensional
overlaps whose derivatives are -1, 0, +1 for a centre and 0, +/-1/2 for a size.
d(chi1inv)/d(fill) comes from varying fill through eff_chi1inv_matrix itself,
via a new fill_override argument, rather than re-deriving Kottke's algebra
here;  is linear in fill by construction so that variation is pure local
algebra with no quadrature in the loop.

The result is step-independent to the bit across a 200x range of the internal
fill step, which is the behaviour an analytic derivative should have and the
old code never had.

The loop now visits only pixels with 0 < fill < 1, which is the discrete form
of a shape derivative being a surface integral.

Magnitudes are still wrong, by a factor that varies with resolution (2.56 at
res 20 against 0.54 at res 40) and with index contrast (2.56 at n=2.5 against
0.74 at n=1.5), so something in the contraction is still structurally wrong
rather than merely scaled. Guarding that get_front_object returned the object
being differentiated, and matching its shiftby in the analytic overlap, changes
nothing -- both were already right.
---
 src/meepgeom.cpp | 387 ++++++++++++++++++++++++++---------------------
 src/meepgeom.hpp |  17 ++-
 2 files changed, 231 insertions(+), 173 deletions(-)

diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 23973a191..44539087b 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -1066,7 +1066,7 @@ void geom_epsilon::eff_chi1inv_row(meep::component c, double chi1inv_row[3], con
 
 void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_matrix,
                                       const meep::volume &v, double tol, int maxeval,
-                                      bool &fallback) {
+                                      bool &fallback, double fill_override) {
   const geometric_object *o;
   material_type mat, mat_behind;
   symm_matrix meps;
@@ -1107,7 +1107,8 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_ma
   pixel.low = vector3_minus(pixel.low, shiftby);
   pixel.high = vector3_minus(pixel.high, shiftby);
 
-  double fill = box_overlap_with_object(pixel, *o, tol, maxeval);
+  double fill =
+      fill_override >= 0.0 ? fill_override : box_overlap_with_object(pixel, *o, tol, maxeval);
 
   material_epsmu(meep::type(c), mat, &meps, chi1inv_matrix);
   symm_matrix eps2, epsinv2;
@@ -1198,6 +1199,29 @@ void geom_epsilon::eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_ma
   sym_matrix_invert(chi1inv_matrix, &meps);
 }
 
+/* The filling fraction this pixel would smooth with. Returns false when the
+   pixel does not straddle an interface between two distinct materials, in
+   which case the shape derivative vanishes there and the caller can skip it --
+   which is what restricts the gradient loop to the boundary shell. */
+bool geom_epsilon::interface_fill(const meep::volume &v, double tol, int maxeval, double &fill,
+                                  const geometric_object **which, vector3 *shift) {
+  const geometric_object *o;
+  material_type mat, mat_behind;
+  vector3 p, shiftby;
+  fill = -1.0;
+
+  if (!get_front_object(v, geometry_tree, p, &o, shiftby, mat, mat_behind)) return false;
+  if (material_type_equal(mat, mat_behind)) return false;
+
+  geom_box pixel = gv2box(v);
+  pixel.low = vector3_minus(pixel.low, shiftby);
+  pixel.high = vector3_minus(pixel.high, shiftby);
+  fill = box_overlap_with_object(pixel, *o, tol, maxeval);
+  if (which) *which = o;
+  if (shift) *shift = shiftby;
+  return true;
+}
+
 static int eps_ever_negative = 0;
 static meep::field_type func_ft = meep::E_stuff;
 
@@ -2945,24 +2969,37 @@ static std::complex<meep::realnum> forward_dft_value(const meep::dft_chunk *ch,
 /* Gradients with respect to a geometric object's centre and size.      */
 /* ------------------------------------------------------------------ */
 
-/* The design gradient finite-differences eff_chi1inv_row with respect to a
-   design weight. eff_chi1inv_row is a function of the *geometry*, so
-   perturbing an object's centre or size instead gives dA/d(parameter) from the
-   same machinery, with no additional timestepping.
+/* The smoothed permittivity depends on the geometry only through two things:
+   the filling fraction of each pixel and the interface normal. So
+
+       d(chi1inv)/d(parameter) = d(chi1inv)/d(fill) * d(fill)/d(parameter)
 
-   One structural difference. A design weight is local, so the design gradient
-   perturbs it inside the point loop. A centre or a size is global, so the
-   perturbation hoists out: the geometry is moved once, every point is visited,
-   and the geometry is restored. That is twelve passes for six parameters, not
-   twelve perturbations per point.
+   with the second factor analytic for an axis-aligned block, since the block
+   is an intersection of three slabs and the overlap volume factorizes:
 
-   The perturbation cannot simply mutate the object, because geom_box_tree
-   caches bounding boxes (`geom_box b, b1, b2` per node, plus one per object).
-   Stale bounds would misreport which object owns a point near the moved
-   boundary -- exactly where a shape derivative lives -- so the tree is rebuilt
-   around each perturbation. */
+       |pixel ∩ block| = prod_i overlap_i(c_i, s_i)
+       overlap_i = clamp(min(c_i + s_i/2, hi_i) - max(c_i - s_i/2, lo_i), 0, dx)
+
+   d(overlap_i)/d(c_i) is -1, 0 or +1 depending on which face lies inside the
+   pixel, and d(overlap_i)/d(s_i) is 0 or +/-1/2. Everything else is a product
+   of the other axes' overlaps.
+
+   The first factor comes from varying `fill` through eff_chi1inv_matrix rather
+   than re-deriving Kottke's algebra here, which keeps this consistent with
+   whatever meep actually does. That variation is pure local algebra: `delta`
+   is linear in fill by construction, and only the final delta -> chi1inv map
+   is nonlinear.
+
+   Perturbing the geometry and re-differencing instead -- which is what this
+   routine used to do -- means differencing box_overlap_with_object, an
+   adaptive quadrature, so its tolerance is amplified by 1/step. That produced
+   a gradient with no stable step regime at all.
+
+   Only pixels straddling the boundary contribute, since d(fill)/d(parameter)
+   vanishes wherever a pixel is wholly inside or wholly outside. That is the
+   discrete form of a shape derivative being a surface integral, and it is why
+   the loop below skips interior pixels. */
 
-// Which scalar of an object is being perturbed.
 enum geom_param {
   GEOM_CENTER_X = 0,
   GEOM_CENTER_Y,
@@ -2972,58 +3009,62 @@ enum geom_param {
   GEOM_SIZE_Z
 };
 
-static double *geom_param_slot(geometric_object *o, int which) {
-  switch (which) {
-    case GEOM_CENTER_X: return &o->center.x;
-    case GEOM_CENTER_Y: return &o->center.y;
-    case GEOM_CENTER_Z: return &o->center.z;
-    default: break;
-  }
-  if (o->which_subclass != geometric_object::BLOCK) return NULL;
-  switch (which) {
-    case GEOM_SIZE_X: return &o->subclass.block_data->size.x;
-    case GEOM_SIZE_Y: return &o->subclass.block_data->size.y;
-    case GEOM_SIZE_Z: return &o->subclass.block_data->size.z;
-    default: return NULL;
-  }
-}
-
-/* Rebuild the box tree after the geometry underneath it has moved. */
-static void geom_rebuild_tree(geom_epsilon *geps, const meep::grid_volume &gv) {
-  if (geps->restricted_tree && geps->restricted_tree != geps->geometry_tree)
-    destroy_geom_box_tree(geps->restricted_tree);
-  destroy_geom_box_tree(geps->geometry_tree);
-  geom_fix_object_list(geps->geometry);
-  geom_box box = gv2box(gv.surroundings());
-  geps->geometry_tree = create_geom_box_tree0(geps->geometry, box);
-  geps->restricted_tree = geps->geometry_tree;
-}
-
-/* One point's contribution: the smoothed chi1inv row at `r`, contracted with
-   the adjoint and forward fields. */
-static std::complex<double> geom_chi1inv_term(geom_epsilon *geps, const meep::grid_volume &gv,
-                                              meep::component adjoint_c, int dir_idx,
-                                              const meep::vec &r, double scale,
-                                              std::complex<meep::realnum> adj,
-                                              std::complex<meep::realnum> fwd) {
-  meep::volume voxel(r);
-  LOOP_OVER_DIRECTIONS(gv.dim, d) {
-    voxel.set_direction_min(d, r.in_direction(d) - 0.5 * gv.inva);
-    voxel.set_direction_max(d, r.in_direction(d) + 0.5 * gv.inva);
-  }
-  double row[3];
-  geps->eff_chi1inv_row(adjoint_c, row, voxel, geps->tol, geps->maxeval);
-  return scale * row[dir_idx] * std::complex<double>(double(adj.real()), double(adj.imag())) *
-         std::complex<double>(double(fwd.real()), double(fwd.imag()));
-}
-
-/* Contract the adjoint and forward fields against the *current* geometry,
-   accumulating one number per frequency. Called twice per parameter, on either
-   side of the perturbation. */
-static void geom_contract(std::complex<double> *out, size_t nf,
-                          std::vector<meep::dft_chunk *> *adjoint_chunks,
-                          std::vector<meep::dft_chunk *> *forward_chunks, double scalegrad,
-                          meep::grid_volume &gv, geom_epsilon *geps) {
+/* Overlap of [lo, hi] with the block's extent along one axis, and its
+   derivatives with respect to that axis's centre and half-size. */
+static void axis_overlap(double lo, double hi, double c, double s, double &overlap,
+                         double &d_dcenter, double &d_dsize) {
+  const double blo = c - 0.5 * s, bhi = c + 0.5 * s;
+  const double left = std::max(lo, blo), right = std::min(hi, bhi);
+  overlap = right - left;
+  if (overlap <= 0) {
+    overlap = 0;
+    d_dcenter = d_dsize = 0;
+    return;
+  }
+  /* Moving the centre moves both faces together; growing the size moves them
+     apart by half each. A face only contributes where it lies strictly inside
+     the pixel -- outside, the clamp pins the overlap and the derivative is
+     zero. This is what makes the derivative one-sided exactly on a pixel
+     edge. */
+  const double lower_inside = (blo > lo) ? 1.0 : 0.0;
+  const double upper_inside = (bhi < hi) ? 1.0 : 0.0;
+  d_dcenter = upper_inside - lower_inside;
+  d_dsize = 0.5 * (upper_inside + lower_inside);
+}
+
+void geometry_addgradient(double *v, size_t nparams, size_t nf,
+                          std::vector<meep::dft_fields *> fields_a,
+                          std::vector<meep::dft_fields *> fields_f, double *frequencies,
+                          double scalegrad, meep::grid_volume &gv, geom_epsilon *geps,
+                          int object_index, int *params, double du) {
+  (void)frequencies;
+  (void)du;
+  if (object_index < 0 || object_index >= geps->geometry.num_items)
+    meep::abort("geometry_addgradient: object index %d out of range (%d objects)", object_index,
+                geps->geometry.num_items);
+  geometric_object *obj = &geps->geometry.items[object_index];
+  if (obj->which_subclass != geometric_object::BLOCK)
+    meep::abort("geometry_addgradient: only blocks are supported");
+
+  const vector3 bc = obj->center;
+  const vector3 bs = obj->subclass.block_data->size;
+  const double centers[3] = {bc.x, bc.y, bc.z};
+  const double sizes[3] = {bs.x, bs.y, bs.z};
+
+  std::vector<meep::dft_chunk *> adjoint_chunks[3], forward_chunks[3];
+  for (int i = 0; i < 3; i++) {
+    for (meep::dft_chunk *c = fields_a[i]->chunks; c; c = c->next_in_dft)
+      adjoint_chunks[i].push_back(c);
+    for (meep::dft_chunk *c = fields_f[i]->chunks; c; c = c->next_in_dft)
+      forward_chunks[i].push_back(c);
+  }
+
+  std::vector<double> local(nf * nparams, 0.0);
+  /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step
+     is well conditioned -- unlike a step in a length, which has to be compared
+     against the pixel. */
+  const double dfill = 1e-6;
+
   for (size_t f_i = 0; f_i < nf; f_i++) {
     for (int ci_adjoint = 0; ci_adjoint < 3; ci_adjoint++) {
       int num_chunks = adjoint_chunks[ci_adjoint].size();
@@ -3039,7 +3080,7 @@ static void geom_contract(std::complex<double> *out, size_t nf,
           if ((num_f == 0) || ((size_t)cur >= num_f)) continue;
           meep::dft_chunk *fwd_chunk = forward_chunks[ci_forward][cur];
           meep::component forward_c = fwd_chunk->c;
-          meep::grid_volume gv_fwd = gv.subvolume(fwd_chunk->is, fwd_chunk->ie, forward_c);
+          if (forward_c != adjoint_c) continue; // diagonal; see note below
 
           int dir_idx;
           switch (meep::component_direction(forward_c)) {
@@ -3051,119 +3092,121 @@ static void geom_contract(std::complex<double> *out, size_t nf,
             default: continue;
           }
 
-          LOOP_OVER_IVECS(gv_adj, adj_chunk->is_old, adj_chunk->ie_old, idx_adj) {
-            IVEC_LOOP_ILOC(gv_adj, ip);
+          meep::ivec loop_is = adj_chunk->persist ? adj_chunk->is_old : adj_chunk->is;
+          meep::ivec loop_ie = adj_chunk->persist ? adj_chunk->ie_old : adj_chunk->ie;
+
+          LOOP_OVER_IVECS(gv_adj, loop_is, loop_ie, idx_adj) {
             IVEC_LOOP_LOC(gv_adj, p);
             std::complex<meep::realnum> adj = adj_chunk->dft[nf * idx_adj + f_i];
-            if (adj == 0.0) continue;
+            std::complex<meep::realnum> fwd = fwd_chunk->dft[nf * idx_adj + f_i];
+            if (adj == 0.0 || fwd == 0.0) continue;
 
-            if (forward_c == adjoint_c) {
-              std::complex<meep::realnum> fwd = fwd_chunk->dft[nf * idx_adj + f_i];
-              double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p.r() : 1;
-              out[f_i] += geom_chi1inv_term(geps, gv, adjoint_c, dir_idx, p, scalegrad * cyl_scale,
-                                            adj, fwd);
+            meep::volume voxel(p);
+            LOOP_OVER_DIRECTIONS(gv.dim, d) {
+              voxel.set_direction_min(d, p.in_direction(d) - 0.5 * gv.inva);
+              voxel.set_direction_max(d, p.in_direction(d) + 0.5 * gv.inva);
             }
-            else {
-              /* Subpixel smoothing is second order because it builds an
-                 effective *tensor*: harmonic averaging for the field component
-                 normal to the interface, arithmetic for the tangential ones.
-                 The off-diagonal entries that produces are non-zero wherever
-                 the interface normal is not along an axis -- which for a shape
-                 derivative is the whole of the signal, since the derivative
-                 lives entirely at the boundary.
-
-                 Contracting them needs the forward field of a different
-                 component, which sits at a different Yee location, so it is
-                 restricted to the two epsilon nodes between the pair and
-                 interpolated there. This mirrors material_grids_addgradient. */
-              meep::ivec fwd_p = ip + gv.iyee_shift(forward_c) - gv.iyee_shift(adjoint_c);
-              meep::ivec unit_a = unit_ivec(gv.dim, component_direction(adjoint_c));
-              meep::ivec unit_f = unit_ivec(gv.dim, component_direction(forward_c));
-              meep::ivec fwd_pa = fwd_p + unit_a * 2;
-              meep::ivec fwd_pf = fwd_p - unit_f * 2;
-              meep::ivec fwd_paf = fwd_p + unit_a * 2 - unit_f * 2;
-
-              meep::ivec fwd_pl[2] = {fwd_p, fwd_pa};
-              meep::ivec fwd_pr[2] = {fwd_pf, fwd_paf};
-              meep::ivec ieps[2] = {(fwd_p + fwd_pf) / 2, (fwd_pa + fwd_paf) / 2};
-
-              for (int node = 0; node < 2; node++) {
-                ptrdiff_t i1 = gv_fwd.index(forward_c, fwd_pl[node]);
-                ptrdiff_t i2 = gv_fwd.index(forward_c, fwd_pr[node]);
-                std::complex<meep::realnum> fwd1 =
-                    ((i1 >= fwd_chunk->N) || (i1 < 0)) ? 0 : fwd_chunk->dft[nf * i1 + f_i];
-                std::complex<meep::realnum> fwd2 =
-                    ((i2 >= fwd_chunk->N) || (i2 < 0)) ? 0 : fwd_chunk->dft[nf * i2 + f_i];
-                std::complex<meep::realnum> fwd_avg =
-                    std::complex<meep::realnum>(0.5, 0) * (fwd1 + fwd2);
-                meep::vec eps1 = gv[ieps[node]];
-                double cyl_scale = (gv.dim == meep::Dcyl) ? eps1.r() : 1;
-                out[f_i] +=
-                    geom_chi1inv_term(geps, gv, adjoint_c, dir_idx, eps1, scalegrad * cyl_scale,
-                                      std::complex<meep::realnum>(0.5, 0) * adj, fwd_avg);
+
+            /* Skip pixels that do not straddle the boundary: d(fill)/d(param)
+               is zero there, so they contribute nothing. */
+            /* `fill` is the fraction of the pixel inside whichever object
+               get_front_object selected. If that is not the object being
+               differentiated, d(fill)/d(our parameters) is not what this pixel
+               responds to -- and taking it anyway gets the sign wrong wherever
+               the front object is the background instead. */
+            double fill;
+            const geometric_object *front = NULL;
+            vector3 shiftby = {0, 0, 0};
+            if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby))
+              continue;
+            if (front != obj) continue;
+            if (fill <= 0.0 || fill >= 1.0) continue;
+
+            /* d(fill)/d(parameter), analytic and separable. */
+            double overlap[3], d_dc[3], d_ds[3], extent[3];
+            bool degenerate = false;
+            for (int ax = 0; ax < 3; ax++) {
+              meep::direction dd = (ax == 0) ? meep::X : (ax == 1) ? meep::Y : meep::Z;
+              const bool resolved = (gv.dim == meep::D3) || (gv.dim == meep::D2 && ax < 2) ||
+                                    (gv.dim == meep::D1 && ax == 2);
+              if (!resolved) {
+                /* A dimension the simulation does not resolve. The block is
+                   effectively infinite along it, so it contributes a factor of
+                   one to the overlap and nothing to the derivative -- note the
+                   block's own size along such an axis is typically zero, so
+                   using it here would make every pixel degenerate. */
+                extent[ax] = 1.0;
+                overlap[ax] = 1.0;
+                d_dc[ax] = d_ds[ax] = 0.0;
+                continue;
               }
+              const double sh = (ax == 0) ? shiftby.x : (ax == 1) ? shiftby.y : shiftby.z;
+              const double lo = voxel.in_direction_min(dd) - sh;
+              const double hi = voxel.in_direction_max(dd) - sh;
+              extent[ax] = hi - lo;
+              axis_overlap(lo, hi, centers[ax], sizes[ax], overlap[ax], d_dc[ax], d_ds[ax]);
+              if (overlap[ax] <= 0) degenerate = true;
+            }
+            if (degenerate) continue;
+
+            double pixel_volume = 1.0;
+            for (int ax = 0; ax < 3; ax++)
+              pixel_volume *= extent[ax];
+
+            /* d(chi1inv)/d(fill), from meep's own tensor assembly. */
+            bool fb_lo = false, fb_hi = false;
+            symm_matrix m_lo, m_hi;
+            geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo,
+                                     std::max(0.0, fill - dfill));
+            geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi,
+                                     std::min(1.0, fill + dfill));
+            if (fb_lo || fb_hi) continue;
+
+            const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02};
+            const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02};
+            const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12};
+            const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12};
+            const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22};
+            const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22};
+            int row_of = 0;
+            switch (meep::component_direction(adjoint_c)) {
+              case meep::X:
+              case meep::R: row_of = 0; break;
+              case meep::Y:
+              case meep::P: row_of = 1; break;
+              case meep::Z: row_of = 2; break;
+              default: continue;
+            }
+            const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2;
+            const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2;
+            const double actual_dfill = std::min(1.0, fill + dfill) - std::max(0.0, fill - dfill);
+            if (actual_dfill <= 0) continue;
+            const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill;
+
+            const std::complex<double> pair =
+                std::complex<double>(double(adj.real()), double(adj.imag())) *
+                std::complex<double>(double(fwd.real()), double(fwd.imag()));
+            const double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p.r() : 1;
+
+            for (size_t ip = 0; ip < nparams; ip++) {
+              const int which = params[ip];
+              const int ax = which % 3;
+              double dfill_dp = (which < 3) ? d_dc[ax] : d_ds[ax];
+              if (dfill_dp == 0.0) continue;
+              for (int other = 0; other < 3; other++)
+                if (other != ax) dfill_dp *= overlap[other];
+              dfill_dp /= pixel_volume;
+
+              /* the leading minus matches get_material_gradient's convention,
+                 which returns -(d row/d parameter) */
+              local[nparams * f_i + ip] -=
+                  scalegrad * cyl_scale * dchi_dfill * dfill_dp * std::real(pair);
             }
           }
         }
       }
     }
   }
-}
-
-/* dJ/d(parameter) for one geometric object.
-
-   `v` receives nf * nparams doubles, frequency-major. `params` names which
-   scalars to differentiate, using the geom_param enumeration.
-
-   Both the diagonal and the off-diagonal terms of the smoothed tensor are
-   contracted; see geom_contract for why the second matters here more than it
-   does for a density gradient. */
-void geometry_addgradient(double *v, size_t nparams, size_t nf,
-                          std::vector<meep::dft_fields *> fields_a,
-                          std::vector<meep::dft_fields *> fields_f, double *frequencies,
-                          double scalegrad, meep::grid_volume &gv, geom_epsilon *geps,
-                          int object_index, int *params, double du) {
-  (void)frequencies;
-  if (object_index < 0 || object_index >= geps->geometry.num_items)
-    meep::abort("geometry_addgradient: object index %d out of range (%d objects)", object_index,
-                geps->geometry.num_items);
-  geometric_object *obj = &geps->geometry.items[object_index];
-
-  std::vector<meep::dft_chunk *> adjoint_chunks[3], forward_chunks[3];
-  for (int i = 0; i < 3; i++) {
-    for (meep::dft_chunk *c = fields_a[i]->chunks; c; c = c->next_in_dft)
-      adjoint_chunks[i].push_back(c);
-    for (meep::dft_chunk *c = fields_f[i]->chunks; c; c = c->next_in_dft)
-      forward_chunks[i].push_back(c);
-  }
-
-  std::vector<std::complex<double> > plus(nf), minus(nf);
-  std::vector<double> local(nf * nparams, 0.0);
-
-  for (size_t ip = 0; ip < nparams; ip++) {
-    double *slot = geom_param_slot(obj, params[ip]);
-    if (!slot)
-      meep::abort("geometry_addgradient: parameter %d is not available on this object", params[ip]);
-    double original = *slot;
-
-    std::fill(minus.begin(), minus.end(), std::complex<double>(0, 0));
-    std::fill(plus.begin(), plus.end(), std::complex<double>(0, 0));
-
-    *slot = original - du;
-    geom_rebuild_tree(geps, gv);
-    geom_contract(minus.data(), nf, adjoint_chunks, forward_chunks, scalegrad, gv, geps);
-
-    *slot = original + du;
-    geom_rebuild_tree(geps, gv);
-    geom_contract(plus.data(), nf, adjoint_chunks, forward_chunks, scalegrad, gv, geps);
-
-    *slot = original;
-    geom_rebuild_tree(geps, gv);
-
-    // (row_1 - row_2) / 2du, matching get_material_gradient's sign convention
-    for (size_t f_i = 0; f_i < nf; f_i++)
-      local[nparams * f_i + ip] = std::real(minus[f_i] - plus[f_i]) / (2 * du);
-  }
 
   meep::sum_to_all(local.data(), v, int(nf * nparams));
 }
diff --git a/src/meepgeom.hpp b/src/meepgeom.hpp
index ce60e1577..dee643ec8 100644
--- a/src/meepgeom.hpp
+++ b/src/meepgeom.hpp
@@ -212,8 +212,23 @@ class geom_epsilon : public meep::material_function {
   virtual bool is_thread_safe() const { return !has_user_materials; }
   bool has_user_materials;
 
+  /* `fill_override >= 0` substitutes that filling fraction for the one this
+     routine would compute from the geometry. The shape derivative uses it: the
+     smoothed tensor depends on the geometry only through `fill` and the
+     interface normal, so d(chi1inv)/d(parameter) factors into
+     d(chi1inv)/d(fill) -- obtained by varying `fill` here, which is pure local
+     algebra with no quadrature -- times d(fill)/d(parameter), which is
+     analytic for an axis-aligned block. Differencing the geometry directly
+     instead means differencing box_overlap_with_object, an adaptive
+     quadrature, which amplifies its tolerance by 1/step. */
   void eff_chi1inv_matrix(meep::component c, symm_matrix *chi1inv_matrix, const meep::volume &v,
-                          double tol, int maxeval, bool &fallback);
+                          double tol, int maxeval, bool &fallback, double fill_override = -1.0);
+
+  /* The filling fraction and interface normal this pixel would smooth with,
+     or fill < 0 if the pixel does not straddle an interface (so the shape
+     derivative is zero there and the pixel can be skipped). */
+  bool interface_fill(const meep::volume &v, double tol, int maxeval, double &fill,
+                      const geometric_object **which = NULL, vector3 *shift = NULL);
 
   void fallback_chi1inv_row(meep::component c, double chi1inv_row[3], const meep::volume &v,
                             double tol, int maxeval);

From f89f8380d7d33f2b543d60e608ba43f99ef342e0 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 11:02:07 -0700
Subject: [PATCH 06/12] wip: fix two support bugs in the analytic geometry
 gradient

Found with a sum rule that needs no fields and no timestepping: summing
d(fill)/d(parameter) times the pixel volume over all pixels must reproduce the
change in the object's own volume -- zero for a translation, and the product of
the other extents for a size.

Two bugs, both of which silently dropped contributions rather than erring.

Strict inequalities in axis_overlap. A face lying exactly on a pixel boundary
satisfied neither `blo > lo` nor `bhi < hi`, so both neighbouring pixels
reported no face inside them and the derivative for that axis vanished
entirely. Round geometry on a round grid hits this constantly: a block of width
1.0 at resolution 20 has faces exactly 10 pixels from its centre. Now half-open,
which assigns such a face to exactly one pixel -- nothing dropped, nothing
double counted, and the derivative one-sided there, which is the truth.

The support test was the wrong set. I argued the support of the shape
derivative is where d(fill)/d(parameter) is non-zero, then implemented
`0 < fill < 1`. Those differ exactly where a face is grid-aligned: no pixel
straddles it, fill is 0 or 1 everywhere along it, and meep applies no smoothing
there at all -- so the pixels carrying that face's derivative were all
excluded. The sum rule showed it plainly: 40 boundary pixels found, which is
exactly the two y faces, with none from the two x faces.

The sum rule itself belongs in the test suite. It is exact, costs no
simulation, and would have caught both of these immediately.
---
 src/meepgeom.cpp | 266 ++++++++++++++++++++++++++++-------------------
 1 file changed, 159 insertions(+), 107 deletions(-)

diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 44539087b..1b39db980 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3022,12 +3022,21 @@ static void axis_overlap(double lo, double hi, double c, double s, double &overl
     return;
   }
   /* Moving the centre moves both faces together; growing the size moves them
-     apart by half each. A face only contributes where it lies strictly inside
-     the pixel -- outside, the clamp pins the overlap and the derivative is
-     zero. This is what makes the derivative one-sided exactly on a pixel
-     edge. */
-  const double lower_inside = (blo > lo) ? 1.0 : 0.0;
-  const double upper_inside = (bhi < hi) ? 1.0 : 0.0;
+     apart by half each. A face contributes only in the pixel that contains it.
+
+     The comparisons are half-open on purpose. With strict inequalities on both
+     sides, a face lying exactly on a pixel boundary belongs to neither
+     neighbour and its contribution vanishes -- silently, and for every pixel
+     along that face, which zeroes the whole derivative for that axis. Round
+     geometry on a round grid hits this constantly: a block of width 1.0 at
+     resolution 20 has faces exactly 10 pixels from its centre.
+
+     Half-open assigns such a face to exactly one of the two pixels, so nothing
+     is dropped and nothing is double counted. The derivative is then one-sided
+     there, which is the truth: the smoothed permittivity has a kink at a pixel
+     boundary and no two-sided derivative exists. */
+  const double lower_inside = (blo >= lo && blo < hi) ? 1.0 : 0.0;
+  const double upper_inside = (bhi > lo && bhi <= hi) ? 1.0 : 0.0;
   d_dcenter = upper_inside - lower_inside;
   d_dsize = 0.5 * (upper_inside + lower_inside);
 }
@@ -3080,7 +3089,7 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
           if ((num_f == 0) || ((size_t)cur >= num_f)) continue;
           meep::dft_chunk *fwd_chunk = forward_chunks[ci_forward][cur];
           meep::component forward_c = fwd_chunk->c;
-          if (forward_c != adjoint_c) continue; // diagonal; see note below
+          meep::grid_volume gv_fwd = gv.subvolume(fwd_chunk->is, fwd_chunk->ie, forward_c);
 
           int dir_idx;
           switch (meep::component_direction(forward_c)) {
@@ -3096,112 +3105,155 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
           meep::ivec loop_ie = adj_chunk->persist ? adj_chunk->ie_old : adj_chunk->ie;
 
           LOOP_OVER_IVECS(gv_adj, loop_is, loop_ie, idx_adj) {
+            IVEC_LOOP_ILOC(gv_adj, ip);
             IVEC_LOOP_LOC(gv_adj, p);
             std::complex<meep::realnum> adj = adj_chunk->dft[nf * idx_adj + f_i];
-            std::complex<meep::realnum> fwd = fwd_chunk->dft[nf * idx_adj + f_i];
-            if (adj == 0.0 || fwd == 0.0) continue;
-
-            meep::volume voxel(p);
-            LOOP_OVER_DIRECTIONS(gv.dim, d) {
-              voxel.set_direction_min(d, p.in_direction(d) - 0.5 * gv.inva);
-              voxel.set_direction_max(d, p.in_direction(d) + 0.5 * gv.inva);
+            if (adj == 0.0) continue;
+
+            /* The smoothed tensor is diagonal in Cartesian axes only where the
+               interface normal lies along an axis, which for a block is true on
+               a face and false on an edge or corner. Those pixels are O(N) of
+               the O(N^2) on the boundary, so their share falls with resolution
+               -- which is the signature of the residual error here. Contracting
+               them needs the forward field of another component, restricted to
+               the two epsilon nodes between the pair. */
+            std::complex<meep::realnum> fwd;
+            meep::vec eps_at = p;
+            double node_weight = 1.0;
+            int num_nodes = 1;
+            meep::vec node_pos[2];
+            std::complex<meep::realnum> node_fwd[2];
+            if (forward_c == adjoint_c) {
+              node_pos[0] = p;
+              node_fwd[0] = fwd_chunk->dft[nf * idx_adj + f_i];
+            }
+            else {
+              num_nodes = 2;
+              node_weight = 0.5;
+              meep::ivec fwd_p = ip + gv.iyee_shift(forward_c) - gv.iyee_shift(adjoint_c);
+              meep::ivec unit_a = unit_ivec(gv.dim, component_direction(adjoint_c));
+              meep::ivec unit_f = unit_ivec(gv.dim, component_direction(forward_c));
+              meep::ivec pl[2] = {fwd_p, fwd_p + unit_a * 2};
+              meep::ivec pr[2] = {fwd_p - unit_f * 2, fwd_p + unit_a * 2 - unit_f * 2};
+              for (int nd = 0; nd < 2; nd++) {
+                ptrdiff_t i1 = gv_fwd.index(forward_c, pl[nd]);
+                ptrdiff_t i2 = gv_fwd.index(forward_c, pr[nd]);
+                std::complex<meep::realnum> f1 =
+                    ((i1 >= fwd_chunk->N) || (i1 < 0)) ? 0 : fwd_chunk->dft[nf * i1 + f_i];
+                std::complex<meep::realnum> f2 =
+                    ((i2 >= fwd_chunk->N) || (i2 < 0)) ? 0 : fwd_chunk->dft[nf * i2 + f_i];
+                node_fwd[nd] = std::complex<meep::realnum>(0.5, 0) * (f1 + f2);
+                node_pos[nd] = gv[(pl[nd] + pr[nd]) / 2];
+              }
             }
 
-            /* Skip pixels that do not straddle the boundary: d(fill)/d(param)
-               is zero there, so they contribute nothing. */
-            /* `fill` is the fraction of the pixel inside whichever object
-               get_front_object selected. If that is not the object being
-               differentiated, d(fill)/d(our parameters) is not what this pixel
-               responds to -- and taking it anyway gets the sign wrong wherever
-               the front object is the background instead. */
-            double fill;
-            const geometric_object *front = NULL;
-            vector3 shiftby = {0, 0, 0};
-            if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby))
-              continue;
-            if (front != obj) continue;
-            if (fill <= 0.0 || fill >= 1.0) continue;
-
-            /* d(fill)/d(parameter), analytic and separable. */
-            double overlap[3], d_dc[3], d_ds[3], extent[3];
-            bool degenerate = false;
-            for (int ax = 0; ax < 3; ax++) {
-              meep::direction dd = (ax == 0) ? meep::X : (ax == 1) ? meep::Y : meep::Z;
-              const bool resolved = (gv.dim == meep::D3) || (gv.dim == meep::D2 && ax < 2) ||
-                                    (gv.dim == meep::D1 && ax == 2);
-              if (!resolved) {
-                /* A dimension the simulation does not resolve. The block is
-                   effectively infinite along it, so it contributes a factor of
-                   one to the overlap and nothing to the derivative -- note the
-                   block's own size along such an axis is typically zero, so
-                   using it here would make every pixel degenerate. */
-                extent[ax] = 1.0;
-                overlap[ax] = 1.0;
-                d_dc[ax] = d_ds[ax] = 0.0;
+            for (int node = 0; node < num_nodes; node++) {
+              fwd = node_fwd[node];
+              if (fwd == 0.0) continue;
+              const meep::vec p_node = node_pos[node];
+
+              meep::volume voxel(p_node);
+              LOOP_OVER_DIRECTIONS(gv.dim, d) {
+                voxel.set_direction_min(d, p_node.in_direction(d) - 0.5 * gv.inva);
+                voxel.set_direction_max(d, p_node.in_direction(d) + 0.5 * gv.inva);
+              }
+
+              /* Skip pixels that do not straddle the boundary: d(fill)/d(param)
+                 is zero there, so they contribute nothing. */
+              /* `fill` is the fraction of the pixel inside whichever object
+                 get_front_object selected. If that is not the object being
+                 differentiated, d(fill)/d(our parameters) is not what this pixel
+                 responds to -- and taking it anyway gets the sign wrong wherever
+                 the front object is the background instead. */
+              double fill;
+              const geometric_object *front = NULL;
+              vector3 shiftby = {0, 0, 0};
+              if (!geps->interface_fill(voxel, geps->tol, geps->maxeval, fill, &front, &shiftby))
                 continue;
+              if (front != obj) continue;
+              if (fill <= 0.0 || fill >= 1.0) continue;
+
+              /* d(fill)/d(parameter), analytic and separable. */
+              double overlap[3], d_dc[3], d_ds[3], extent[3];
+              bool degenerate = false;
+              for (int ax = 0; ax < 3; ax++) {
+                meep::direction dd = (ax == 0) ? meep::X : (ax == 1) ? meep::Y : meep::Z;
+                const bool resolved = (gv.dim == meep::D3) || (gv.dim == meep::D2 && ax < 2) ||
+                                      (gv.dim == meep::D1 && ax == 2);
+                if (!resolved) {
+                  /* A dimension the simulation does not resolve. The block is
+                     effectively infinite along it, so it contributes a factor of
+                     one to the overlap and nothing to the derivative -- note the
+                     block's own size along such an axis is typically zero, so
+                     using it here would make every pixel degenerate. */
+                  extent[ax] = 1.0;
+                  overlap[ax] = 1.0;
+                  d_dc[ax] = d_ds[ax] = 0.0;
+                  continue;
+                }
+                const double sh = (ax == 0) ? shiftby.x : (ax == 1) ? shiftby.y : shiftby.z;
+                const double lo = voxel.in_direction_min(dd) - sh;
+                const double hi = voxel.in_direction_max(dd) - sh;
+                extent[ax] = hi - lo;
+                axis_overlap(lo, hi, centers[ax], sizes[ax], overlap[ax], d_dc[ax], d_ds[ax]);
+                if (overlap[ax] <= 0) degenerate = true;
               }
-              const double sh = (ax == 0) ? shiftby.x : (ax == 1) ? shiftby.y : shiftby.z;
-              const double lo = voxel.in_direction_min(dd) - sh;
-              const double hi = voxel.in_direction_max(dd) - sh;
-              extent[ax] = hi - lo;
-              axis_overlap(lo, hi, centers[ax], sizes[ax], overlap[ax], d_dc[ax], d_ds[ax]);
-              if (overlap[ax] <= 0) degenerate = true;
-            }
-            if (degenerate) continue;
-
-            double pixel_volume = 1.0;
-            for (int ax = 0; ax < 3; ax++)
-              pixel_volume *= extent[ax];
-
-            /* d(chi1inv)/d(fill), from meep's own tensor assembly. */
-            bool fb_lo = false, fb_hi = false;
-            symm_matrix m_lo, m_hi;
-            geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo,
-                                     std::max(0.0, fill - dfill));
-            geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi,
-                                     std::min(1.0, fill + dfill));
-            if (fb_lo || fb_hi) continue;
-
-            const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02};
-            const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02};
-            const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12};
-            const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12};
-            const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22};
-            const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22};
-            int row_of = 0;
-            switch (meep::component_direction(adjoint_c)) {
-              case meep::X:
-              case meep::R: row_of = 0; break;
-              case meep::Y:
-              case meep::P: row_of = 1; break;
-              case meep::Z: row_of = 2; break;
-              default: continue;
-            }
-            const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2;
-            const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2;
-            const double actual_dfill = std::min(1.0, fill + dfill) - std::max(0.0, fill - dfill);
-            if (actual_dfill <= 0) continue;
-            const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill;
-
-            const std::complex<double> pair =
-                std::complex<double>(double(adj.real()), double(adj.imag())) *
-                std::complex<double>(double(fwd.real()), double(fwd.imag()));
-            const double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p.r() : 1;
-
-            for (size_t ip = 0; ip < nparams; ip++) {
-              const int which = params[ip];
-              const int ax = which % 3;
-              double dfill_dp = (which < 3) ? d_dc[ax] : d_ds[ax];
-              if (dfill_dp == 0.0) continue;
-              for (int other = 0; other < 3; other++)
-                if (other != ax) dfill_dp *= overlap[other];
-              dfill_dp /= pixel_volume;
-
-              /* the leading minus matches get_material_gradient's convention,
-                 which returns -(d row/d parameter) */
-              local[nparams * f_i + ip] -=
-                  scalegrad * cyl_scale * dchi_dfill * dfill_dp * std::real(pair);
-            }
+              if (degenerate) continue;
+
+              double pixel_volume = 1.0;
+              for (int ax = 0; ax < 3; ax++)
+                pixel_volume *= extent[ax];
+
+              /* d(chi1inv)/d(fill), from meep's own tensor assembly. */
+              bool fb_lo = false, fb_hi = false;
+              symm_matrix m_lo, m_hi;
+              geps->eff_chi1inv_matrix(adjoint_c, &m_lo, voxel, geps->tol, geps->maxeval, fb_lo,
+                                       std::max(0.0, fill - dfill));
+              geps->eff_chi1inv_matrix(adjoint_c, &m_hi, voxel, geps->tol, geps->maxeval, fb_hi,
+                                       std::min(1.0, fill + dfill));
+              if (fb_lo || fb_hi) continue;
+
+              const double lo_row[3] = {m_lo.m00, m_lo.m01, m_lo.m02};
+              const double hi_row[3] = {m_hi.m00, m_hi.m01, m_hi.m02};
+              const double lo_row1[3] = {m_lo.m01, m_lo.m11, m_lo.m12};
+              const double hi_row1[3] = {m_hi.m01, m_hi.m11, m_hi.m12};
+              const double lo_row2[3] = {m_lo.m02, m_lo.m12, m_lo.m22};
+              const double hi_row2[3] = {m_hi.m02, m_hi.m12, m_hi.m22};
+              int row_of = 0;
+              switch (meep::component_direction(adjoint_c)) {
+                case meep::X:
+                case meep::R: row_of = 0; break;
+                case meep::Y:
+                case meep::P: row_of = 1; break;
+                case meep::Z: row_of = 2; break;
+                default: continue;
+              }
+              const double *lo_r = (row_of == 0) ? lo_row : (row_of == 1) ? lo_row1 : lo_row2;
+              const double *hi_r = (row_of == 0) ? hi_row : (row_of == 1) ? hi_row1 : hi_row2;
+              const double actual_dfill = std::min(1.0, fill + dfill) - std::max(0.0, fill - dfill);
+              if (actual_dfill <= 0) continue;
+              const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill;
+
+              const std::complex<double> pair =
+                  std::complex<double>(double(adj.real()), double(adj.imag())) *
+                  std::complex<double>(double(fwd.real()), double(fwd.imag()));
+              const double cyl_scale = (gv.dim == meep::Dcyl) ? 2 * p_node.r() : 1;
+
+              for (size_t ip = 0; ip < nparams; ip++) {
+                const int which = params[ip];
+                const int ax = which % 3;
+                double dfill_dp = (which < 3) ? d_dc[ax] : d_ds[ax];
+                if (dfill_dp == 0.0) continue;
+                for (int other = 0; other < 3; other++)
+                  if (other != ax) dfill_dp *= overlap[other];
+                dfill_dp /= pixel_volume;
+
+                /* the leading minus matches get_material_gradient's convention,
+                   which returns -(d row/d parameter) */
+                local[nparams * f_i + ip] -=
+                    node_weight * scalegrad * cyl_scale * dchi_dfill * dfill_dp * std::real(pair);
+              }
+            } // node
           }
         }
       }

From 7467099d9457145e7cd09ec862830e71277f6800 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 11:20:05 -0700
Subject: [PATCH 07/12] wip: isolate the defect to d(fill)/d(parameter), with a
 field-free check

MEEP_CHECK_DFILL compares the analytic d(fill)/d(parameter) against a direct
measurement, per pixel, with no fields and no chi1inv in the path. The
measurement calls box_overlap_with_object on the object itself, so it needs no
geometry tree and cannot be tripped by the stale bounding boxes that moving an
object leaves behind.

Validating the diagnostic before trusting it: its reference scale comes out at
exactly 0.5/Delta -- 10.0 at resolution 20 and 20.0 at resolution 40 -- which
is what a size parameter must give when one face sweeps a pixel. Two earlier
versions of this same check were wrong and I did trust them. One used a
quarter-pixel step, large enough that fill saturated and the reference merely
measured the clamp. The other read fill back through interface_fill after
moving the object without rebuilding the tree, so get_front_object returned
garbage.

With a sound reference, the analytic derivative is off by up to 50% on some
pixels. So the defect is in d(fill)/d(parameter) specifically, not in
d(chi1inv)/d(fill), not in the field contraction, and not in the off-diagonal
terms -- those are exactly zero for an axis-aligned block anyway, since
normal_to_fixed_object returns the nearest face's normal and the Kottke tensor
stays diagonal in Cartesian axes.

Next: dump the per-pixel comparison and find which pixels disagree. The error
is a constant 5.0 in absolute terms at both resolutions, which does not fit a
corner-pixel explanation and needs the pattern to identify.
---
 src/meepgeom.cpp | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 1b39db980..3dad3ed1c 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3069,6 +3069,8 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
   }
 
   std::vector<double> local(nf * nparams, 0.0);
+  double fill_worst = 0.0, fill_scale = 0.0;
+  size_t fill_checked = 0;
   /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step
      is well conditioned -- unlike a step in a length, which has to be compared
      against the pixel. */
@@ -3248,6 +3250,36 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
                   if (other != ax) dfill_dp *= overlap[other];
                 dfill_dp /= pixel_volume;
 
+                if (f_i == 0 && node == 0 && getenv("MEEP_CHECK_DFILL")) {
+                  /* Measure d(fill)/dp directly. box_overlap_with_object takes
+                     the object, so it needs no geometry tree and cannot be
+                     tripped by the stale bounding boxes that moving an object
+                     leaves behind. */
+                  const double hh = 0.01 / gv.a;
+                  double *slot = (which < 3) ? ((ax == 0)   ? &obj->center.x
+                                                : (ax == 1) ? &obj->center.y
+                                                            : &obj->center.z)
+                                             : ((ax == 0)   ? &obj->subclass.block_data->size.x
+                                                : (ax == 1) ? &obj->subclass.block_data->size.y
+                                                            : &obj->subclass.block_data->size.z);
+                  const double orig = *slot;
+                  const geom_box pixbox = gv2box(voxel);
+                  *slot = orig + hh;
+                  geom_fix_object_list(geps->geometry);
+                  const double f_hi =
+                      box_overlap_with_object(pixbox, *obj, geps->tol, geps->maxeval);
+                  *slot = orig - hh;
+                  geom_fix_object_list(geps->geometry);
+                  const double f_lo =
+                      box_overlap_with_object(pixbox, *obj, geps->tol, geps->maxeval);
+                  *slot = orig;
+                  geom_fix_object_list(geps->geometry);
+                  const double numeric = (f_hi - f_lo) / (2 * hh);
+                  fill_checked++;
+                  fill_scale = std::max(fill_scale, std::abs(numeric));
+                  fill_worst = std::max(fill_worst, std::abs(numeric - dfill_dp));
+                }
+
                 /* the leading minus matches get_material_gradient's convention,
                    which returns -(d row/d parameter) */
                 local[nparams * f_i + ip] -=
@@ -3260,6 +3292,9 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
     }
   }
 
+  if (getenv("MEEP_CHECK_DFILL"))
+    master_printf("DFILL checked=%zu worst_abs_err=%.6e scale=%.6e\n", fill_checked, fill_worst,
+                  fill_scale);
   meep::sum_to_all(local.data(), v, int(nf * nparams));
 }
 

From f5745a6bb530b38d7d94c7a36ae765451e423e88 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 11:49:29 -0700
Subject: [PATCH 08/12] geometry gradient: d(fill)/d(parameter) verified exact
 to 1e-12

The field-free per-pixel check now agrees with a direct measurement to 1.2e-12
against a scale of 10, at resolutions 20 and 40. So the analytic fill
derivative is right, and whatever remains is elsewhere.

Getting there needed one correction and one piece of structure.

The voxel-edge warning was inverted. Pixel centres are at integer multiples of
dx, so pixel edges are at half-integers: a face at a pixel *centre* is
straddled and smooths normally, and a face at a pixel *edge* is straddled by
nothing. The warning fired on integers -- the good case -- and the half-pixel
offset I had introduced to "fix" an earlier sweep moved faces from centres onto
edges. I had been measuring the degenerate configuration and calling it well
conditioned.

The structure is that this is component-dependent. Yee components sit half a
pixel apart, so one component's grid puts a given face mid-pixel while
another's puts it exactly on an edge, simultaneously. There is no placement
that is well conditioned for every component at once. Faces have to avoid both
integer and half-integer multiples of dx, so the safe offset is a quarter
pixel, not a half. With that, every component straddles every face and the
derivative is two-sided everywhere.

Where a face does land on some component's pixel edge, the disagreement is not
an error in either quantity: the overlap is clamped at zero or full, so the
analytic value is the correct one-sided derivative while a central difference
reports the average of the two sides. Both are right; they answer different
questions.

Still wrong: the assembled gradient, by a factor that varies (0.77, 1.15, 2.18
across resolution and index). With d(fill)/d(parameter) exact, that isolates to
d(chi1inv)/d(fill) or the field contraction. The same field-free treatment
applies to the former and is the next step.
---
 python/adjoint/geometry_gradient.py | 15 +++++++++++----
 src/meepgeom.cpp                    | 10 ++++++++++
 2 files changed, 21 insertions(+), 4 deletions(-)

diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
index 0176676b8..4e332714e 100644
--- a/python/adjoint/geometry_gradient.py
+++ b/python/adjoint/geometry_gradient.py
@@ -147,17 +147,24 @@ def check_faces_off_voxel_edges(sim: mp.Simulation, obj) -> None:
         if size[axis] == 0:
             continue
         for face in (center[axis] - size[axis] / 2, center[axis] + size[axis] / 2):
+            # Pixel centres sit at integer multiples of dx, so pixel *edges*
+            # are at half-integers. A face at a pixel centre is straddled by
+            # that pixel and smooths normally; a face at a pixel edge is
+            # straddled by nothing, every neighbouring pixel is full or empty,
+            # and the derivative there is one-sided.
             offset = abs((face / dx) - round(face / dx))
-            if offset < 1e-6:
+            if abs(offset - 0.5) < 1e-6:
                 warnings.warn(
                     f"The {letter} face of "
                     f"{getattr(obj, 'name', None) or 'this object'} at "
-                    f"{face:.6g} lies on a voxel edge, where the smoothed "
+                    f"{face:.6g} lies on a pixel edge, where the smoothed "
                     "permittivity has a kink and the derivative with respect "
                     "to position is one-sided. The reported gradient will "
                     "depend on the finite-difference step rather than "
-                    f"converging. Offset the geometry by about {dx / 2:.6g} "
-                    "along " + letter + ", or change the resolution.",
+                    f"converging. Shift it by about {dx / 2:.6g} along "
+                    + letter
+                    + " so the face lands on a pixel centre, or change the "
+                    "resolution.",
                     RuntimeWarning,
                     stacklevel=3,
                 )
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 3dad3ed1c..b5e9d35cd 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3070,6 +3070,7 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
 
   std::vector<double> local(nf * nparams, 0.0);
   double fill_worst = 0.0, fill_scale = 0.0;
+  int fill_printed = 0;
   size_t fill_checked = 0;
   /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step
      is well conditioned -- unlike a step in a length, which has to be compared
@@ -3278,6 +3279,15 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
                   fill_checked++;
                   fill_scale = std::max(fill_scale, std::abs(numeric));
                   fill_worst = std::max(fill_worst, std::abs(numeric - dfill_dp));
+                  if (std::abs(numeric - dfill_dp) > 1e-9 && fill_printed < 12) {
+                    fill_printed++;
+                    master_printf("  BAD p%d at (%.4f,%.4f) fill=%.4f ov=(%.4f,%.4f,%.4f) "
+                                  "dc=(%.2f,%.2f,%.2f) ds=(%.2f,%.2f,%.2f) vol=%.5f "
+                                  "analytic=%.4f numeric=%.4f\n",
+                                  which, p_node.x(), p_node.y(), fill, overlap[0], overlap[1],
+                                  overlap[2], d_dc[0], d_dc[1], d_dc[2], d_ds[0], d_ds[1], d_ds[2],
+                                  pixel_volume, dfill_dp, numeric);
+                  }
                 }
 
                 /* the leading minus matches get_material_gradient's convention,

From bfeb5cf7e6a90b1e9031555ac2cb3a4a0b9587d3 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 12:41:45 -0700
Subject: [PATCH 09/12] geometry shape derivative works: 0.1-0.9% against
 finite differences

Verified at two resolutions and two index contrasts, with the residual being
the finite difference's own truncation -- a step sweep shows it still
converging toward the adjoint at h = 0.02 pixels.

The last several rounds of apparent failure were a broken test, not broken
code. t_ratio.py perturbed the block's centre to form the finite difference
while reading the adjoint's *size* entry, left over from an earlier edit that
switched it to size for a sum-rule check. It had been comparing
d/d(center.y) against d/d(size.y). A test that compares centre to centre
agrees to 0.19%.

What actually needed fixing, in order of discovery:

  - the geometry finite difference, which differenced an adaptive quadrature
    and so amplified its tolerance by 1/step, replaced by the analytic chain
    rule d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp;
  - strict inequalities in axis_overlap, which zeroed a whole axis whenever a
    face lay on a pixel boundary;
  - the support test, which used 0 < fill < 1 rather than d(fill)/dp != 0;
  - an inverted voxel-edge warning: pixel centres are the well conditioned
    case and pixel edges the degenerate one, and Yee components sit half a
    pixel apart, so faces must avoid both integer and half-integer multiples
    of dx -- a quarter-pixel offset is safe for every component at once.

Both operator factors are independently verified rather than inferred:
d(fill)/dp agrees with a direct measurement to 1.2e-12, and d(chi1inv)/d(fill)
matches the closed form implied by Kottke's algebra -- -(eps1-eps2)/eps_avg^2
for the tangential component -- to six digits at two different fill values.
---
 src/meepgeom.cpp | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index b5e9d35cd..ee88f0cf6 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3071,6 +3071,7 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
   std::vector<double> local(nf * nparams, 0.0);
   double fill_worst = 0.0, fill_scale = 0.0;
   int fill_printed = 0;
+  int dchi_printed = 0;
   size_t fill_checked = 0;
   /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step
      is well conditioned -- unlike a step in a length, which has to be compared
@@ -3237,6 +3238,18 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
               if (actual_dfill <= 0) continue;
               const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill;
 
+            if (f_i == 0 && node == 0 && getenv("MEEP_CHECK_DCHI") && dchi_printed < 14) {
+              /* Kottke assembles the normal component harmonically, so
+                 chi1inv_nn = fill/eps1 + (1-fill)/eps2 is linear in fill and
+                 its derivative is exactly 1/eps1 - 1/eps2, independent of
+                 fill. The tangential component is arithmetic in eps, giving
+                 -(eps1-eps2)/eps_avg^2. Both are known in closed form, so this
+                 validates the instrument rather than assuming it. */
+              dchi_printed++;
+              master_printf("  DCHI c=%s dir=%d fill=%.4f dchi_dfill=%.6f\n",
+                            meep::component_name(adjoint_c), dir_idx, fill, dchi_dfill);
+            }
+
               const std::complex<double> pair =
                   std::complex<double>(double(adj.real()), double(adj.imag())) *
                   std::complex<double>(double(fwd.real()), double(fwd.imag()));

From 4a62d4c341513908a10126f70b0ba8923c319ae8 Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 28 Aug 2026 13:15:49 -0700
Subject: [PATCH 10/12] adjoint: gradients with respect to a geometric object's
 centre and size

An object opts in with `differentiable=['center', 'size']` and its gradient
comes back under its `name`, alongside the design gradient. Costs no extra
simulation.

Subpixel smoothing makes the permittivity depend on geometry only through each
pixel's filling fraction and the interface normal, so

    d(chi1inv)/dp = d(chi1inv)/d(fill) * d(fill)/dp

The first factor is obtained by varying `fill` through eff_chi1inv_matrix
itself, via a new fill_override argument, rather than re-deriving Kottke's
algebra -- `delta` is linear in fill by construction, so that variation is pure
local algebra. The second is analytic: a block is an intersection of slabs, so
the pixel overlap factorizes into clamped one-dimensional overlaps.

Differencing the geometry instead would mean differencing
box_overlap_with_object, an adaptive quadrature, which amplifies its tolerance
by 1/step and has no stable step regime at all.

Only pixels the boundary passes through contribute, since d(fill)/dp vanishes
wherever a pixel is wholly inside or outside. That is the discrete form of a
shape derivative being a surface integral.

Two conventions the caller has to know, both documented and one enforced.
Subpixel smoothing is required, and the gradient refuses without it rather than
returning a number, because an unsmoothed permittivity is a step function of
position. And object faces want to sit a quarter pixel clear of pixel edges:
Yee components are half a pixel apart, so a face at a pixel centre for one
component is on a pixel edge for another, where nothing straddles it and the
derivative is one-sided. A warning fires when that is detected.

Verified against finite differences at 0.1-0.9% across resolution and index
contrast, with the residual being the finite difference's own truncation -- a
step sweep shows it still converging toward the adjoint. Both operator factors
are also checked independently: d(fill)/dp against a direct measurement to
1.2e-12, and d(chi1inv)/d(fill) against the closed form Kottke's algebra
implies, to six digits.
---
 NEWS.md                                     |  11 +
 doc/docs/Python_Tutorials/Adjoint_Solver.md |  66 ++++++
 python/Makefile.am                          |   2 +
 python/tests/test_geometry_gradient.py      | 234 ++++++++++++++++++++
 src/meepgeom.cpp                            |  58 -----
 5 files changed, 313 insertions(+), 58 deletions(-)
 create mode 100644 python/tests/test_geometry_gradient.py

diff --git a/NEWS.md b/NEWS.md
index 0bcde63bd..21d36a170 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -2,6 +2,17 @@
 
 ## Meep 1.35.0 (in progress)
 
+* Adjoint solver: geometric objects can be differentiated with respect to their
+  `center` and `size`, so an optimizer can be asked whether a reflector should
+  move or a spacer should lengthen. An object opts in with
+  `differentiable=['center', 'size']` and its gradient appears under its `name`.
+  It costs no extra simulation. Only pixels the object's boundary passes through
+  contribute, which is the discrete form of a shape derivative being a surface
+  integral. Requires subpixel smoothing, and refuses to run without it. Object
+  faces should be kept a quarter pixel clear of pixel edges, since Yee
+  components sit half a pixel apart and a face on some component's pixel edge
+  has only a one-sided derivative; a warning is issued when this is detected.
+
 * Adjoint solver: sources can now be differentiated alongside the design
   regions. A source opts in by naming parameters, as in
   `differentiable=['beam_w0', 'beam_x0']`, and its gradient appears in the
diff --git a/doc/docs/Python_Tutorials/Adjoint_Solver.md b/doc/docs/Python_Tutorials/Adjoint_Solver.md
index b5245afe1..86c47f87a 100644
--- a/doc/docs/Python_Tutorials/Adjoint_Solver.md
+++ b/doc/docs/Python_Tutorials/Adjoint_Solver.md
@@ -492,6 +492,72 @@ the run was doubled.
 Finally, a source inside or near a PML is rejected: its adjoint field is absorbed,
 so the gradient would come back finite, smooth, and wrong.
 
+Differentiating With Respect To Geometry
+----------------------------------------
+
+The design region is not the only thing that can move. An ordinary geometric
+object can be differentiated too — where a reflector sits, how long a spacer is
+— by naming the parameters on the object itself:
+
+```py
+reflector = mp.Block(
+    center=mp.Vector3(0, 2.0),
+    size=mp.Vector3(4, 0.3),
+    material=silicon,
+    differentiable=["center", "size"],
+    name="reflector",
+)
+
+value, grad = opt([rho])
+grad["design"]                # as before
+grad["reflector"]["center"]   # (3,), or (nfreq, 3) for several frequencies
+grad["reflector"]["size"]
+```
+
+This is the same pattern as the design weights and the Gaussian beam
+parameters: finite-difference a cheap analytic map and contract it against the
+adjoint field, at no cost in extra timestepping. Subpixel smoothing makes the
+permittivity depend on the geometry only through each pixel's filling fraction
+and the interface normal, so
+
+$$\frac{\partial \chi^{-1}}{\partial p} = \frac{\partial \chi^{-1}}{\partial f}\,\frac{\partial f}{\partial p}$$
+
+with $\partial f/\partial p$ analytic — a block is an intersection of slabs, so
+the pixel overlap factorizes — and only pixels the boundary passes through
+contributing anything. That is the discrete form of a shape derivative being a
+surface integral.
+
+`'center'` and `'size'` are accepted here and rejected on a *source*. That is
+not an inconsistency: a source's cotangent is gathered over a fixed set of grid
+points, so moving it changes which points it occupies rather than the
+amplitudes applied to them, whereas the permittivity is a function of position
+and moving an object is exactly what a derivative with respect to position
+means.
+
+### Two things to get right
+
+**Subpixel smoothing must be on.** It is on by default (`eps_averaging=True`),
+and the gradient refuses to run without it rather than returning a number.
+Without smoothing the permittivity is a step function of position — nothing
+changes until a boundary crosses a pixel edge, then it changes by the full
+material contrast — and a finite difference of that is not a derivative.
+
+**Keep object faces off the pixel edges.** Pixel centres lie at integer
+multiples of the pixel, so pixel edges lie at half-integers, and Yee components
+sit half a pixel apart from one another. A face at a pixel centre for one field
+component therefore lies exactly on a pixel edge for another, where no pixel
+straddles it and the derivative is one-sided. A **quarter-pixel** offset clears
+both, so every component straddles every face:
+
+```py
+offset = 0.25 / resolution
+block = mp.Block(center=mp.Vector3(0, y0 + offset), ...)
+```
+
+Round geometry on a round grid lands on edges constantly — a block 1.0 wide at
+resolution 20 has faces exactly 10 pixels from its centre — so this is worth
+doing deliberately. A warning is issued when a face is detected on an edge.
+
 Broadband Waveguide Mode Converter with Minimum Feature Size
 ------------------------------------------------------------
 
diff --git a/python/Makefile.am b/python/Makefile.am
index 2ba73a947..efe5e4ff3 100644
--- a/python/Makefile.am
+++ b/python/Makefile.am
@@ -37,6 +37,7 @@ ADJOINT_TESTS =                                \
     $(TEST_DIR)/test_adjoint_symmetry.py       \
     $(TEST_DIR)/test_adjoint_protocol.py       \
     $(TEST_DIR)/test_angular_spectrum.py       \
+    $(TEST_DIR)/test_geometry_gradient.py      \
     $(TEST_DIR)/test_source_gradient.py        \
     $(TEST_DIR)/test_adjoint_jax.py
 
@@ -101,6 +102,7 @@ TESTS =                                        \
     $(TEST_DIR)/test_simulation.py             \
     $(TEST_DIR)/test_special_kz.py             \
     $(TEST_DIR)/test_source.py                 \
+    $(TEST_DIR)/test_geometry_gradient.py      \
     $(TEST_DIR)/test_source_gradient.py        \
     $(TEST_DIR)/test_stop_when_flux_decayed.py \
     $(TEST_DIR)/test_subpixel_3d.py            \
diff --git a/python/tests/test_geometry_gradient.py b/python/tests/test_geometry_gradient.py
new file mode 100644
index 000000000..3cc142d8e
--- /dev/null
+++ b/python/tests/test_geometry_gradient.py
@@ -0,0 +1,234 @@
+"""Tests for adjoint gradients with respect to a geometric object's shape.
+
+Two conventions run through all of these, and both cost real debugging time to
+establish.
+
+Object faces are placed at *quarter*-pixel offsets. Yee components sit half a
+pixel apart, so a face at a pixel centre for one component lies exactly on a
+pixel edge for another. On an edge no pixel straddles the face: the filling
+fraction is 0 or 1 on both sides, meep applies no smoothing there, and the
+derivative is one-sided. A quarter-pixel offset avoids both integer and
+half-integer multiples of the pixel, so every component straddles every face
+and the derivative is two-sided everywhere.
+
+Run lengths are pinned rather than left to `stop_when_dft_decayed`, for the
+same reason as the source-gradient tests: an adaptive stop makes a perturbed
+run end at a different time, and the difference scales with the perturbation.
+"""
+
+import unittest
+import warnings
+
+import numpy as np
+from autograd import numpy as npa
+
+import meep as mp
+import meep.adjoint as mpa
+
+RES = 20
+FCEN = 1.0
+RUN = 200.0
+CELL = mp.Vector3(8, 6)
+SRC_C = mp.Vector3(-2.0, 0)
+MON_C = mp.Vector3(2.0, 0)
+RHO = 0.5 * np.ones(64)
+FD_STEP = 1e-3
+
+
+def quarter_pixel(resolution: float) -> float:
+    """The offset that keeps a face clear of every component's pixel edges."""
+    return 0.25 / resolution
+
+
+def _design_region(sim):
+    """An inert design region; OptimizationProblem requires one."""
+    grid = mp.MaterialGrid(
+        mp.Vector3(8, 8), mp.air, mp.Medium(index=1.5), grid_type="U_MEAN"
+    )
+    region = mpa.DesignRegion(
+        grid, volume=mp.Volume(center=mp.Vector3(0, -2.0), size=mp.Vector3(0.4, 0.4))
+    )
+    sim.geometry = list(sim.geometry) + [
+        mp.Block(center=region.center, size=region.size, material=grid)
+    ]
+    return region
+
+
+def _problem(block, res=RES, component=mp.Hz, eps_averaging=True):
+    sim = mp.Simulation(
+        cell_size=CELL,
+        resolution=res,
+        boundary_layers=[mp.PML(1.0)],
+        sources=[
+            mp.Source(
+                mp.GaussianSource(FCEN, fwidth=0.2), component=component, center=SRC_C
+            )
+        ],
+        geometry=[block],
+        eps_averaging=eps_averaging,
+        force_complex_fields=True,
+    )
+    region = _design_region(sim)
+    monitor = mpa.FourierFields(
+        sim, mp.Volume(center=MON_C, size=mp.Vector3(0, 0)), component
+    )
+    return mpa.OptimizationProblem(
+        simulation=sim,
+        objective_functions=[lambda f: npa.sum(npa.abs(f) ** 2)],
+        objective_arguments=[monitor],
+        design_regions=[region],
+        frequencies=[FCEN],
+        minimum_run_time=RUN,
+        maximum_run_time=RUN,
+    )
+
+
+def _block(center, size, index=2.5, differentiable=None, res=RES):
+    """A block whose faces avoid every component's pixel edges."""
+    offset = quarter_pixel(res)
+    return mp.Block(
+        center=mp.Vector3(center.x, center.y + offset),
+        size=mp.Vector3(size.x + 2 * offset, size.y),
+        material=mp.Medium(index=index),
+        differentiable=differentiable,
+        name="scatterer",
+    )
+
+
+class TestDifferentiableFlag(unittest.TestCase):
+    """Validation of `differentiable=` on a geometric object."""
+
+    def test_accepts_center_and_size(self):
+        block = mp.Block(
+            center=mp.Vector3(),
+            size=mp.Vector3(1, 1),
+            differentiable=["center", "size"],
+            name="b",
+        )
+        self.assertEqual(block.differentiable, ("center", "size"))
+        self.assertEqual(block.name, "b")
+
+    def test_default_is_not_differentiable(self):
+        self.assertEqual(mp.Block(size=mp.Vector3(1, 1)).differentiable, ())
+
+    def test_rejects_unknown_parameter(self):
+        with self.assertRaisesRegex(ValueError, "not a differentiable parameter"):
+            mp.Block(size=mp.Vector3(1, 1), differentiable=["waist"])
+
+    def test_rejects_bare_string(self):
+        with self.assertRaisesRegex(ValueError, "list of parameter names"):
+            mp.Block(size=mp.Vector3(1, 1), differentiable="center")
+
+    def test_center_and_size_are_allowed_here_unlike_for_sources(self):
+        # A source's cotangent is gathered over a fixed set of grid points, so
+        # moving it changes which points it occupies rather than the amplitudes
+        # applied to them, and 'center' is rejected there. Epsilon is a
+        # function of position, so for geometry the same name is exactly what a
+        # derivative with respect to position means.
+        with self.assertRaises(ValueError):
+            mp.Source(
+                mp.GaussianSource(FCEN, fwidth=0.2),
+                component=mp.Ez,
+                center=SRC_C,
+                differentiable=["center"],
+            )
+        mp.Block(size=mp.Vector3(1, 1), differentiable=["center"])
+
+
+class TestGuards(unittest.TestCase):
+    def test_refuses_without_subpixel_smoothing(self):
+        # Without smoothing the permittivity is a step function of position:
+        # nothing changes until a boundary crosses a pixel edge, then it
+        # changes by the full contrast. A finite difference of that is not a
+        # derivative, and returning a number would hide it.
+        block = _block(mp.Vector3(), mp.Vector3(1.0, 0.8), differentiable=["center"])
+        opt = _problem(block, eps_averaging=False)
+        with self.assertRaisesRegex(ValueError, "subpixel smoothing"):
+            opt([RHO])
+
+    def test_warns_when_a_face_lies_on_a_pixel_edge(self):
+        # Pixel centres are at integer multiples of dx, so edges are at
+        # half-integers. A face there is straddled by nothing.
+        half = 0.5 / RES
+        block = mp.Block(
+            center=mp.Vector3(0, half),
+            size=mp.Vector3(1.0, 0.8),
+            material=mp.Medium(index=2.5),
+            differentiable=["center"],
+            name="s",
+        )
+        opt = _problem(block)
+        with warnings.catch_warnings(record=True) as caught:
+            warnings.simplefilter("always")
+            opt([RHO])
+        self.assertTrue(
+            any("pixel edge" in str(w.message) for w in caught),
+            "expected a warning about a face on a pixel edge",
+        )
+
+
+class TestShapeDerivative(unittest.TestCase):
+    SIZE = mp.Vector3(1.0, 0.8)
+
+    def _value(self, center, size, res=RES, index=2.5):
+        opt = _problem(_block(center, size, index, res=res), res=res)
+        return float(np.asarray(opt([RHO], need_gradient=False)[0]).item())
+
+    def _adjoint(self, names, res=RES, index=2.5):
+        opt = _problem(
+            _block(mp.Vector3(), self.SIZE, index, differentiable=names, res=res),
+            res=res,
+        )
+        _, grad = opt([RHO])
+        return grad["scatterer"]
+
+    def test_center_gradient(self):
+        for label, res, index in (
+            ("baseline", RES, 2.5),
+            ("higher resolution", 40, 2.5),
+            ("lower contrast", RES, 1.5),
+        ):
+            with self.subTest(label):
+                gradient = self._adjoint(["center"], res, index)["center"]
+                adjoint = float(np.real(np.atleast_1d(gradient)[1]))
+                d = mp.Vector3(y=FD_STEP)
+                reference = (
+                    self._value(d, self.SIZE, res, index)
+                    - self._value(mp.Vector3() - d, self.SIZE, res, index)
+                ) / (2 * FD_STEP)
+                self.assertLess(
+                    abs(adjoint - reference) / abs(reference),
+                    2e-2,
+                    f"{label}: adjoint {adjoint} vs finite difference {reference}",
+                )
+
+    def test_size_gradient(self):
+        adjoint = float(np.real(np.atleast_1d(self._adjoint(["size"])["size"])[1]))
+        d = mp.Vector3(y=FD_STEP)
+        reference = (
+            self._value(mp.Vector3(), self.SIZE + d)
+            - self._value(mp.Vector3(), self.SIZE - d)
+        ) / (2 * FD_STEP)
+        self.assertLess(abs(adjoint - reference) / abs(reference), 2e-2)
+
+    def test_translation_through_a_uniform_medium_is_zero(self):
+        # A block of the background material is not there at all, so moving it
+        # cannot change anything. Exact, and independent of finite-difference
+        # error, which makes it sharper than the comparisons above.
+        gradient = self._adjoint(["center"], index=1.0)["center"]
+        self.assertLess(float(np.max(np.abs(np.atleast_1d(gradient)))), 1e-12)
+
+    def test_returns_one_entry_per_named_parameter(self):
+        gradient = self._adjoint(["center", "size"])
+        self.assertEqual(set(gradient), {"center", "size"})
+        for name in ("center", "size"):
+            self.assertEqual(np.shape(np.atleast_1d(gradient[name])), (3,))
+
+    def test_no_flagged_object_leaves_the_return_shape_alone(self):
+        opt = _problem(_block(mp.Vector3(), self.SIZE))
+        _, gradient = opt([RHO])
+        self.assertNotIsInstance(gradient, dict)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index ee88f0cf6..1b39db980 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3069,10 +3069,6 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
   }
 
   std::vector<double> local(nf * nparams, 0.0);
-  double fill_worst = 0.0, fill_scale = 0.0;
-  int fill_printed = 0;
-  int dchi_printed = 0;
-  size_t fill_checked = 0;
   /* `fill` is dimensionless and O(1) and `delta` is linear in it, so this step
      is well conditioned -- unlike a step in a length, which has to be compared
      against the pixel. */
@@ -3238,18 +3234,6 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
               if (actual_dfill <= 0) continue;
               const double dchi_dfill = (hi_r[dir_idx] - lo_r[dir_idx]) / actual_dfill;
 
-            if (f_i == 0 && node == 0 && getenv("MEEP_CHECK_DCHI") && dchi_printed < 14) {
-              /* Kottke assembles the normal component harmonically, so
-                 chi1inv_nn = fill/eps1 + (1-fill)/eps2 is linear in fill and
-                 its derivative is exactly 1/eps1 - 1/eps2, independent of
-                 fill. The tangential component is arithmetic in eps, giving
-                 -(eps1-eps2)/eps_avg^2. Both are known in closed form, so this
-                 validates the instrument rather than assuming it. */
-              dchi_printed++;
-              master_printf("  DCHI c=%s dir=%d fill=%.4f dchi_dfill=%.6f\n",
-                            meep::component_name(adjoint_c), dir_idx, fill, dchi_dfill);
-            }
-
               const std::complex<double> pair =
                   std::complex<double>(double(adj.real()), double(adj.imag())) *
                   std::complex<double>(double(fwd.real()), double(fwd.imag()));
@@ -3264,45 +3248,6 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
                   if (other != ax) dfill_dp *= overlap[other];
                 dfill_dp /= pixel_volume;
 
-                if (f_i == 0 && node == 0 && getenv("MEEP_CHECK_DFILL")) {
-                  /* Measure d(fill)/dp directly. box_overlap_with_object takes
-                     the object, so it needs no geometry tree and cannot be
-                     tripped by the stale bounding boxes that moving an object
-                     leaves behind. */
-                  const double hh = 0.01 / gv.a;
-                  double *slot = (which < 3) ? ((ax == 0)   ? &obj->center.x
-                                                : (ax == 1) ? &obj->center.y
-                                                            : &obj->center.z)
-                                             : ((ax == 0)   ? &obj->subclass.block_data->size.x
-                                                : (ax == 1) ? &obj->subclass.block_data->size.y
-                                                            : &obj->subclass.block_data->size.z);
-                  const double orig = *slot;
-                  const geom_box pixbox = gv2box(voxel);
-                  *slot = orig + hh;
-                  geom_fix_object_list(geps->geometry);
-                  const double f_hi =
-                      box_overlap_with_object(pixbox, *obj, geps->tol, geps->maxeval);
-                  *slot = orig - hh;
-                  geom_fix_object_list(geps->geometry);
-                  const double f_lo =
-                      box_overlap_with_object(pixbox, *obj, geps->tol, geps->maxeval);
-                  *slot = orig;
-                  geom_fix_object_list(geps->geometry);
-                  const double numeric = (f_hi - f_lo) / (2 * hh);
-                  fill_checked++;
-                  fill_scale = std::max(fill_scale, std::abs(numeric));
-                  fill_worst = std::max(fill_worst, std::abs(numeric - dfill_dp));
-                  if (std::abs(numeric - dfill_dp) > 1e-9 && fill_printed < 12) {
-                    fill_printed++;
-                    master_printf("  BAD p%d at (%.4f,%.4f) fill=%.4f ov=(%.4f,%.4f,%.4f) "
-                                  "dc=(%.2f,%.2f,%.2f) ds=(%.2f,%.2f,%.2f) vol=%.5f "
-                                  "analytic=%.4f numeric=%.4f\n",
-                                  which, p_node.x(), p_node.y(), fill, overlap[0], overlap[1],
-                                  overlap[2], d_dc[0], d_dc[1], d_dc[2], d_ds[0], d_ds[1], d_ds[2],
-                                  pixel_volume, dfill_dp, numeric);
-                  }
-                }
-
                 /* the leading minus matches get_material_gradient's convention,
                    which returns -(d row/d parameter) */
                 local[nparams * f_i + ip] -=
@@ -3315,9 +3260,6 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
     }
   }
 
-  if (getenv("MEEP_CHECK_DFILL"))
-    master_printf("DFILL checked=%zu worst_abs_err=%.6e scale=%.6e\n", fill_checked, fill_worst,
-                  fill_scale);
   meep::sum_to_all(local.data(), v, int(nf * nparams));
 }
 

From 10442ae7ad9034f671555e89bbcd27af750d8c3c Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Mon, 31 Aug 2026 13:26:12 -0700
Subject: [PATCH 11/12] Clamp geometry monitors to the cell

`mp.inf` is the idiomatic way to say "spans the cell", and it is how one writes
a layer in a stratified stack.  It is 1e20 rather than a flag, and
`_fit_volume_to_simulation` passes it straight through, so the padded monitor
built around a differentiable object of infinite extent asked for a DFT volume
1e20 wide and meep failed with "impossible(?) looping boundaries".

Clamp each padded extent to the cell before fitting.

Found while building a grating coupler with a metal reflector, where the
reflector and every other layer in the stack are naturally written with an
infinite in-plane extent, so this is hit immediately rather than as a corner
case.
---
 python/adjoint/geometry_gradient.py    | 16 +++++++++++---
 python/tests/test_geometry_gradient.py | 30 ++++++++++++++++++++++++++
 2 files changed, 43 insertions(+), 3 deletions(-)

diff --git a/python/adjoint/geometry_gradient.py b/python/adjoint/geometry_gradient.py
index 4e332714e..d35490a25 100644
--- a/python/adjoint/geometry_gradient.py
+++ b/python/adjoint/geometry_gradient.py
@@ -211,10 +211,20 @@ def install_geometry_monitors(
         # where the field is continuous across the boundary and a great deal
         # where it is not -- the normal E field jumps by the index contrast.
         pad = MONITOR_PAD_PIXELS / sim.resolution
+
+        # `mp.inf` is how one writes "spans the cell", and is the natural way to
+        # describe a layer -- but it is 1e20, not a flag, and
+        # `_fit_volume_to_simulation` passes it straight through. A DFT volume
+        # of that extent fails inside meep with "impossible(?) looping
+        # boundaries", so clamp to the cell before fitting.
+        cell = sim.cell_size
         padded = mp.Vector3(
-            obj.size.x + 2 * pad if obj.size.x else 0.0,
-            obj.size.y + 2 * pad if obj.size.y else 0.0,
-            obj.size.z + 2 * pad if obj.size.z else 0.0,
+            *(
+                min(s + 2 * pad, c) if s else 0.0
+                for s, c in zip(
+                    (obj.size.x, obj.size.y, obj.size.z), (cell.x, cell.y, cell.z)
+                )
+            )
         )
         volume = sim._fit_volume_to_simulation(
             mp.Volume(center=obj.center, size=padded)
diff --git a/python/tests/test_geometry_gradient.py b/python/tests/test_geometry_gradient.py
index 3cc142d8e..3a3d24595 100644
--- a/python/tests/test_geometry_gradient.py
+++ b/python/tests/test_geometry_gradient.py
@@ -224,6 +224,36 @@ def test_returns_one_entry_per_named_parameter(self):
         for name in ("center", "size"):
             self.assertEqual(np.shape(np.atleast_1d(gradient[name])), (3,))
 
+    def test_object_with_an_infinite_extent(self):
+        # `mp.inf` is the idiomatic way to write a layer that spans the cell,
+        # and it is 1e20 rather than a flag. A monitor built around one without
+        # clamping fails inside meep with "impossible(?) looping boundaries",
+        # which is what a stratified stack -- a metal reflector under a grating,
+        # say -- runs into immediately.
+        offset = quarter_pixel(RES)
+
+        def slab(dy, differentiable=None):
+            return mp.Block(
+                center=mp.Vector3(0, dy + offset),
+                size=mp.Vector3(mp.inf, 0.4),
+                material=mp.Medium(index=2.5),
+                differentiable=differentiable,
+                name="slab",
+            )
+
+        opt = _problem(slab(0.0, differentiable=["center"]))
+        _, grad = opt([RHO])
+        adjoint = float(np.real(np.atleast_1d(grad["slab"]["center"])[1]))
+
+        def value(dy):
+            return float(
+                np.asarray(_problem(slab(dy))([RHO], need_gradient=False)[0]).item()
+            )
+
+        reference = (value(FD_STEP) - value(-FD_STEP)) / (2 * FD_STEP)
+        self.assertGreater(abs(reference), 1e-6, "objective must respond to the move")
+        self.assertLess(abs(adjoint - reference) / abs(reference), 2e-2)
+
     def test_no_flagged_object_leaves_the_return_shape_alone(self):
         opt = _problem(_block(mp.Vector3(), self.SIZE))
         _, gradient = opt([RHO])

From 3f8d8ffd6e42a172fbfe4af22c111a67714c607d Mon Sep 17 00:00:00 2001
From: Alec Hammond <alechammond@meta.com>
Date: Fri, 4 Sep 2026 10:31:33 -0700
Subject: [PATCH 12/12] Pair forward and adjoint chunks by position, not by
 list index

geometry_addgradient took the forward chunk as forward_chunks[ci][cur], using
the adjoint chunk's index in its own list.  That assumes the two lists have the
same length and ordering; they need not, which is exactly why
matching_dft_chunk exists a few functions above for material_grids_addgradient.
This routine was added later and did not use it.

The consequence is worse than a mixed-up gradient.  The mispaired chunk is then
indexed with the adjoint chunk's offsets: node_fwd[0] used idx_adj, an offset
into gv_adj, to index fwd_chunk->dft, with no bounds check -- the two lookups
immediately below it already had one.  Pair spatially and route the lookup
through forward_dft_value, which checks the point is inside the chunk and the
index inside its array.
---
 src/meepgeom.cpp | 20 ++++++++++++++++----
 1 file changed, 16 insertions(+), 4 deletions(-)

diff --git a/src/meepgeom.cpp b/src/meepgeom.cpp
index 1b39db980..96c595350 100644
--- a/src/meepgeom.cpp
+++ b/src/meepgeom.cpp
@@ -3085,9 +3085,17 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
         meep::grid_volume gv_adj = gv.subvolume(adj_chunk->is, adj_chunk->ie, adjoint_c);
 
         for (int ci_forward = 0; ci_forward < 3; ci_forward++) {
-          size_t num_f = forward_chunks[ci_forward].size();
-          if ((num_f == 0) || ((size_t)cur >= num_f)) continue;
-          meep::dft_chunk *fwd_chunk = forward_chunks[ci_forward][cur];
+          /* Pair the forward chunk with this adjoint chunk *spatially*.
+             Indexing both lists by `cur` assumes they are the same length and
+             in the same order, and they need not be -- which is why
+             matching_dft_chunk exists a few functions up, for exactly this
+             purpose in material_grids_addgradient.  A foreign chunk here does
+             not merely mix the gradient up: it then gets indexed with this
+             chunk's offsets and read past the end of its array. */
+          if (forward_chunks[ci_forward].empty()) continue;
+          meep::dft_chunk *fwd_chunk =
+              matching_dft_chunk(forward_chunks[ci_forward], adj_chunk);
+          if (!fwd_chunk) continue;
           meep::component forward_c = fwd_chunk->c;
           meep::grid_volume gv_fwd = gv.subvolume(fwd_chunk->is, fwd_chunk->ie, forward_c);
 
@@ -3125,7 +3133,11 @@ void geometry_addgradient(double *v, size_t nparams, size_t nf,
             std::complex<meep::realnum> node_fwd[2];
             if (forward_c == adjoint_c) {
               node_pos[0] = p;
-              node_fwd[0] = fwd_chunk->dft[nf * idx_adj + f_i];
+              /* `idx_adj` is an offset into gv_adj, so it cannot be used to
+                 index the forward chunk.  Look the value up by location, which
+                 bounds-checks it as the two lookups below already did. */
+              node_fwd[0] =
+                  forward_dft_value(fwd_chunk, fwd_chunk->dft, gv_fwd, gv, ip, nf, f_i);
             }
             else {
               num_nodes = 2;