Smooth a MaterialGrid object's own boundary, and make its interior average consistent for graded grids

[jin-castle]

do_averaging=True promises subpixel smoothing for a MaterialGrid, but two
regimes never got it right, and both matter for using the C++ smoothing path
with externally projected designs. One commit per defect, a third that adds
symmetry-invariance and anisotropic-media coverage this path never had, and a
fourth that makes the interface-in-voxel question explicit — and the quadrature
cheaper — now that it can be asked without introducing a discontinuity. The
harmonic-average commit is separable into its own PR on request.

1. The object's own boundary is never smoothed. In 3d the top and bottom
   faces of a design slab — the interfaces that set the vertical confinement of
   every strip-waveguide design — stay staircased at every resolution, no
   matter what do_averaging is set to.
2. The interior average is inconsistent for graded (unprojected, beta=0)
   grids. The harmonic term integrates the inverse of a binary mixture
   rather than of the graded medium the rest of the code assembles, which
   leaves the assembled epsilon discontinuous in the design weights at the
   |∇u| < 1e-8 gate — and the adjoint gradient wrong by 6–47 % for weights
   projected upstream (e.g. SSP), the exact composition where a low-resolution
   design is interpolated to the Yee grid through the meep API.

Commit 1: the material-grid fallback is the only smoothing path that cannot see an object boundary

There are two paths. eff_chi1inv_matrix() calls get_front_object() and, if
that succeeds, averages two materials across the object's boundary using its
normal and the pixel's filling fraction. If it fails, fallback_chi1inv_row()
integrates along a level set instead.

get_front_object() rejected any pixel touching a variable material:

if ((o1 && is_variable(o1->material)) || (o2 && is_variable(o2->material)) || ...)
  return false;

So a material-grid object's own boundary never reached the first path. And the
fallback takes its normal from ∇u alone for a grid (src/meepgeom.cpp:1406)
while a user material function gets it from normal_vector() and integrates
chi1p1() (:1409) — both of which sample the actual material and therefore
already see the object boundary. Only the material-grid branch is blind to
it.

For an (n, n, 1) grid ∇u has no z component at all. The slab's z faces are
not in the design level set; they are the object's boundary. Nothing in either
path was smoothing them. Those faces do not exist in 2d, and every adjoint test
upstream is 2d or cylindrical, which is why this survived.

Why the front material may be a grid and the back material may not

The average needs one material that is constant over the pixel to act as the
background, and one boundary with a normal and a filling fraction. mat_behind
must therefore stay constant — a variable material there is still rejected.
mat_front need not be: the boundary being averaged is the object's own level
set, and the grid can be evaluated at a point before averaging.

Which point matters. The pixel center can sit on the far side of the boundary,
and that is exactly where a MaterialGrid is evaluated on its mirrored
extension. For a planar cut of a pixel of size h with inside fraction fill,
the centroid of the inside part lies (1-fill)*h/2 from the center along the
inward normal. libctl does not fix the normal's orientation, so both candidates
are tried and the one point_in_fixed_objectp() accepts is kept.

Interior pixels must go back to the fallback. A pixel with fill = 1 has no
boundary in it; averaging it would return the unsmoothed material and silently
disable level-set smoothing across the whole design region. Those are handed
back.

Changes (commit 1)

1. Narrow the rejection in get_front_object() to what actually breaks the
   average, and move it after the front/behind assignment so it can distinguish
   the two. A material grid is admissible as the front object's own material
   (not as default_material, which has no boundary to average) when
   do_averaging is set; front_is_matgrid tells the caller to localize it.
2. inside_sample_point() — evaluate the grid at the centroid of the inside
   part of the pixel, per the derivation above.
3. Boundary pixels only. fill >= 1 - 1e-6 returns fallback = true;
   fill <= 1e-6 takes the no-average path.
4. Defer the metal test for a front-side grid until after localization.
   is_metal() reads the grid's medium, which before get_material_pt()
   still holds the previous pixel's interpolation. mat_behind is never a grid,
   so its test stays where it was.

Commit 2: integrate 1/eps of the medium the code actually assembles

matgrid_inveps_func() (and its complex-integration twin) integrated

(1 - u_proj) / eps1 + u_proj / eps2

— the inverse of a binary mixture, valid only where u_proj is 0 or 1. The
medium every other code path assembles is graded:
eps(u) = (1-u)*eps1 + u*eps2. For intermediate u_proj the two differ by the
arithmetic–harmonic gap, and that gap does not vanish as |∇u| → 0. So the
two branches of the fallback disagreed at the gate: below |∇u| < 1e-8 a pixel
gets the isotropic pointwise 1/eps(uval); just above it, the anisotropic
combination (1/meps)·I + (minveps - 1/meps)·nnᵀ with a normal-direction
epsilon that is O(1) different (at u = 0.5, index 1.44 vs 3.48: 7.09
pointwise vs 3.53 anisotropic). The assembled epsilon is discontinuous in the
design weights wherever the design field has a flat spot — which is the norm
for weights projected upstream at beta=inf (SSP) and handed to a beta=0
grid, and happens at every critical point of a smooth design field too. The
adjoint evaluates fields on the closed branch and dA/du on the open one, and
the mismatch shows up directly against a directional finite difference.

Integrating 1/eps(u_proj) instead sends minveps → 1/meps in the
zero-gradient limit, so the anisotropic correction vanishes continuously and
the two branches agree. For projected (binary) u_proj the integrand is
unchanged — projected-grid users see the same numbers, which is what
test_adjoint_solver and test_material_grid passing unchanged confirms.

Two consequences handled in the same commit:

• The voxel model's linear extrapolation u = uval + |∇u|·x is saturated to
  [0, 1]. The interpolated weight itself never leaves that range — only the
  extrapolation does, at a steep transition next to a saturated node — and an
  out-of-range weight extrapolates epsilon through zero, where 1/eps has a
  pole (the old linear-in-u integrand merely went unphysical there; the new one
  went infinite).
• One caveat: a grid whose two media straddle eps = 0 (metal–dielectric) now
  meets a genuine 1/eps singularity inside the integral where the old
  integrand stayed finite. get_front_object() already declines to average
  metals; the fallback inherits the general position that averaging metals is
  not meaningful.

Measured (commit 2)

Adjoint vs directional finite difference along the adjoint gradient; 3d
strip-waveguide design region, EigenmodeCoefficient objective, design 17×17×1,
resolution 12, grid beta=0, eps_averaging=True, do_averaging=True:

design weights	dp	before	after
SSP-projected (beta=inf upstream), pulled inside [0,1]	1e-3	6.6 %	0.094 %
same	3e-3	47.5 %	0.003 %
conic-filtered (smooth)	1e-3	1.8 %	0.033 %
SSP-projected, do_averaging=False (control)	1e-3	0.000 %	0.000 %

The before-column's strong dp-dependence is the discontinuity itself: the
secant straddles jumps. The smooth field's 1.8 % was the same defect at the
field's isolated critical points (∇u ≈ 0), not FD noise.

This is the mechanism behind #2757. The reproduction there is a constant
design region — a ∇u = 0 plateau sitting exactly on the discontinuous gate —
and a structure-level probe of that configuration (no FDTD: perturb the
weights, read back mp.Dielectric) shows the step is gone. Uniform u = 0.9
design, do_averaging=True, perturbation 1e-5:

measurement	before	after
random perturbation, relative move in epsilon	1.97e-01	5.1e-06
smoothed vs unsmoothed epsilon at ∇u = 0	10.81 vs 8.32	identical
the step across beta = 0 / 1 / 4	0.20 / 0.18 / 0.056	~1e-6
the step across resolution 10–60	0.22–0.25 (mesh-independent)	~5e-6

A 1e-5 nudge now moves epsilon by 1e-5, at every beta and every
resolution — which is what made #2757's finite difference irreconcilable with
its adjoint. I'd suggest confirming on the original script there before
closing the issue from this PR.

Commit 4: an explicit interface-distance test, as a quadrature skip

Under the same linear voxel model the quadrature uses, the η level set
crosses a voxel only if |uval − η| < |∇u|·rad. tanh saturates to ±1 in
double precision once |β·(u − η)| exceeds ~19, so when the whole voxel sits
at least that far on one side of η, u_proj is constant across it to
machine precision and the average is the pointwise material — return it
directly instead of integrating a constant.

The condition is written as a saturation width, |uval − η| − |∇u|·rad > 20/β, not as a hard |d| < rad cut: a hard cut at finite β would
reintroduce a (small) step at the threshold — the disease commit 2 just
removed. At β = inf it reduces to |uval − η| > |∇u|·rad exactly, and
β = 0 (a graded medium with no level-set interface) keeps the full average.

Measured: the assembled epsilon is unchanged (bit-identical at β = 8,
≤ 5.3e-15 at β = inf — quadrature roundoff vs. a direct division on the
same constant), and set_epsilon on a 3d random binary design drops
1.41 s → 1.29 s at β = inf — a worst case, since a random binary design
puts an interface next to nearly every voxel; contiguous designs skip far
more.

Tests

TestObjectBoundarySmoothing in python/tests/test_subpixel_3d.py. Filling
weights with a constant removes the in-plane level set entirely, leaving the
block's six faces as the only material interfaces — so whatever smoothing the
grid gets there must be the smoothing an ordinary Medium block of the same
epsilon gets.

• z profiles across the slab's faces agree with the equivalent plain block to
  1e-8;
• the do_averaging on/off difference exceeds 0.5, which guards the first
  assertion against passing because both sides are staircased — before commit 1
  that difference was identically zero;
• a diagonal ramp at beta=inf keeps a genuine level-set interface in the
  interior and must show an on/off difference above 1.0 (measured 3.91) — this
  is the guard that interior pixels are not swallowed by the fill = 1
  boundary path. Diagonal, because for an axis-aligned linear ramp the voxel
  windows uval ± |∇u|·rad tile the u-axis exactly (2·rad is the pixel), so
  at most one column per component straddles the interface and the signal is an
  alignment accident;
• the same ramp at beta=0 is a graded medium with no interface anywhere, and
  must show an on/off difference below 0.1 (measured 0.0024) — before commit 2
  it read 3.9, the binary harmonic mix applied to a medium with no interface in
  any voxel.

A third commit adds coverage for two regimes this path never had tests for:

• symmetries invariance — the smoothed epsilon assembly is identical with
  and without mp.Mirror(mp.Y) for a mirror-symmetric design (measured
  3.6e-15, at beta=inf and beta=0). Structure only: the adjoint gradient
  under symmetries has its own independent, pre-existing defect that nothing
  here touches.
• Anisotropic media — a uniform grid of anisotropic media gets the same
  analytic boundary average as a plain block of the interpolated anisotropic
  epsilon (identical to the last bit — same code path), with an on/off guard
  showing the average engages; and in the interior, where the level-set
  fallback deliberately declines anisotropic pixels (its 1d normal average is
  scalar), do_averaging is pinned as a finite no-op rather than a
  half-applied average.

Full suites on a clean build of this branch: test_subpixel_3d.py 11/11,
test_material_grid.py 3/3 (transmission identical to all printed digits),
test_adjoint_solver.py all pass.

One CI note: the MPI jobs initially tripped on
test_gradient_matches_fd_with_smoothing — an existing test from #3263 that
holds fd/adjoint to 1%. The excess is the chunk-dependence of the adjoint
gradient that #3274 fixes: it is np-dependent (1.3% at np=2, 7.3% at np=3,
serial exact) and vanishes at every rank count on this branch merged with
#3274's head. Until #3274 lands, the parallel assertion is bounded at 0.1 —
still far below the 70–100% disagreement the test guards against — and the
serial assertion keeps its 1%.

Also checked against the pending #3277 (70161188), since both touch
meepgeom.cpp: the merge is conflict-free, and on the combined build both
sides' suites pass — this branch's 11/11 plus #3277's
test_adjoint_symmetric_grids (3/3) and test_adjoint_adjacent_grids (5/5).
Two composition results worth recording: a MaterialGrid overlapped with a
mirrored copy of itself (U_MEAN, the #1984 pattern) assembles the same
smoothed epsilon as a single grid carrying the explicitly averaged weights to
2.7e-15; and #3277's symmetric-grids test rerun with do_averaging=True
passes 3/3 (adjoint/FD = 1.0 ± 2e-3) — a configuration that needs both #3277's
shared-variable attribution and this PR's graded harmonic average to hold.

What would falsify this

If the uniform-grid z profile does not match the equivalent block, the
localization point or the filling fraction is wrong. If the beta=inf ramp's
on/off difference collapses, interior pixels are being averaged against
fill = 1 and level-set smoothing is gone. If the beta=0 ramp's difference
grows back, the anisotropic correction is again surviving without an interface
in the voxel. All three are asserted, so the tests are the falsifier.

Not addressed

#1965 (edge artifacts that vanish when the grid is default_material
rather than a Block) may have the same root as commit 1, but I have not
confirmed it.

The two commits are separable — happy to split the harmonic-integrand change
into its own PR if preferred.

cc @smartalecH @stevengj

[lxvm]

This would be a nice feature to have since the workaround of defining a 3d material grid to smooth out the normal boundaries of a 2d design region with finite thickness could cause you to run out of memory when doing parallel meep calculations since you would have to load a large 3d material grid onto each process

[smartalecH]

This isn't related to the core PR goal (smooth MG boundaries), right? The native cpp MG smoothing code indeed needs to be revamped, but I think we should do that separately.

[smartalecH]

I'm not sure this part will properly get reflected in the $\frac{dA}{du}$ part of the adjoint's recombination step (currently uses a finite difference over the smoothing routine... but not sure this part of the smoothing routine). Can we do some additional gradient checks to see?