[BUG] The natural boundary condition is not working on the internal boundary

[NengLu]

The natural boundary condition is not working on the internal boundary.

Here is a simple test: uniform viscosity, 0 body force, traction boundary condition on internal boundary, no slip on other boundaries.

uw3

import underworld3 as uw
import numpy as np
import sympy
xmin, xmax = 0, 1
ymin, yint,ymax = 0,0.5,1
xres, yres = 32,32
yresa,yresb = 16,16
mesh = mesh = uw.meshing.BoxInternalBoundary(elementRes=(xres,yres),zelementRes=(yresa,yresb),minCoords=(xmin,ymin),maxCoords=(xmax, ymax),zintCoord=yint,degree=1,qdegree=2)
v = uw.discretisation.MeshVariable("V", mesh, mesh.dim, degree=2,continuous=True)
p = uw.discretisation.MeshVariable("P", mesh, 1, degree=1,continuous=True)
stokes = uw.systems.Stokes(mesh, velocityField=v, pressureField=p)
stokes.constitutive_model = uw.constitutive_models.ViscousFlowModel
stokes.bodyforce = sympy.Matrix([0., 0.])
stokes.constitutive_model.Parameters.shear_viscosity_0 = 1.0
stokes.saddle_preconditioner = 1.0 / stokes.constitutive_model.Parameters.shear_viscosity_0
for wall in ("Left", "Right", "Bottom", "Top"):
    stokes.add_essential_bc((0.0, 0.0), wall)
stokes.add_natural_bc(sympy.Matrix([1.0, 0.0]), "Internal")
stokes.tolerance = 1.0e-6
stokes.petsc_options["ksp_rtol"] = 1.0e-6
stokes.petsc_options["ksp_atol"] = 1.0e-6
stokes.petsc_options["snes_converged_reason"] = None
stokes.petsc_options["snes_monitor_short"] = None
stokes.solve(zero_init_guess=False)
print("u just below y=0.5:", uw.function.evaluate(v,np.array([0.5,0.499])))
print("u just above y=0.5:", uw.function.evaluate(v,np.array([0.5,0.501])))

The solver did not converge.

same case using FEniCS

from dolfin import *
import numpy as np
mesh = UnitSquareMesh(32, 32)
OUTER     = 1
INTERFACE = 2
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1, 0)
class OuterBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary
class Interface(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.5)   
OuterBoundary().mark(boundaries, OUTER)
Interface().mark(boundaries, INTERFACE)
dx  = Measure("dx",  domain=mesh)
#ds  = Measure("ds",  domain=mesh, subdomain_data=boundaries)  # exterior facets
dS  = Measure("dS",  domain=mesh, subdomain_data=boundaries)  # interior facets
P2 = VectorElement("P", mesh.ufl_cell(), 2)
P1 = FiniteElement("P", mesh.ufl_cell(), 1)
W  = FunctionSpace(mesh, MixedElement([P2, P1]))
(u, p) = TrialFunctions(W)
(v, q) = TestFunctions(W)
mu = Constant(1.0)
f  = Constant((0.0, 0.0))
def eps(u):
    return sym(grad(u))
a = (inner(2*mu*eps(u), eps(v))*dx
     - p*div(v)*dx
     - q*div(u)*dx )
L = inner(f, v)*dx
# Traction on interior interface y = 0.5 
t = Constant((1.0, 0.0))  # horizontal shear traction
L += dot(t, v('+')) * dS(INTERFACE)
# No-slip on outer boundary
bc_noslip = DirichletBC(W.sub(0), Constant((0.0, 0.0)), boundaries, OUTER)
# remove the null space
bc_pressure = DirichletBC(W.sub(1), Constant(0.0),"x[0] < DOLFIN_EPS && x[1] < DOLFIN_EPS","pointwise")
bcs = [bc_noslip, bc_pressure]
w = Function(W)
solve(a == L, w, bcs,
      solver_parameters={"linear_solver":  "mumps"})
(u_sol, p_sol) = w.split(deepcopy=True)
print("u just below y=0.5:", u_sol(0.5, 0.499))
print("u just above y=0.5:", u_sol(0.5, 0.501))

The output is:

u just below y=0.5: [5.65310750e-02 6.15222229e-09]
u just above y=0.5: [5.65310750e-02 6.15222231e-09]

[lmoresi]

@NengLu — confirmed and root cause located. The bug isn't internal-boundary-specific; it's a label-population issue in `BoxInternalBoundary` that affects every boundary on that mesh (`Top`, `Bottom`, `Left`, `Right`, and `Internal`).

Diagnosis

I reproduced your case and got u ≈ 2.5e-19 everywhere — the SNES "converges" in 0 iterations because the initial zero velocity already satisfies F < 1e-11. That can only mean the natural BC contribution is silently absent from the residual.

To check whether `add_natural_bc` works at all, I ran two control cases:

1. Vanilla `StructuredQuadBox` + natural BC on Top → works correctly. Initial F-norm 0.131, converges to a Couette-like profile, u_top ≈ -0.15. So the natural BC machinery itself is functional.

2. `BoxInternalBoundary` + natural BC on Top (exterior!) → same failure as your `Internal` case — 0 iterations, u ≈ 2.5e-19.

So the issue isn't "natural BC on internal boundary doesn't work". It's "natural BC doesn't work at all on `BoxInternalBoundary`".

Where the difference lies

Comparing label populations on the two mesh classes:

Label	`StructuredQuadBox`	`BoxInternalBoundary`
`Top`	63 points (facets + vertex closure)	32 points (facets only)
`All_Boundaries`	128	0
`Null_Boundary`	1089	0

The vanilla mesh labels not just the height-1 facets but also their bounding vertices (the closure), and populates the consolidated `All_Boundaries` / `Null_Boundary` labels. `BoxInternalBoundary` labels only the facets and leaves the consolidated labels empty.

PETSc's natural-BC residual integration depends on these auxiliary labels in several paths (boundary closure for FE basis evaluation; `All_Boundaries` for some quadrature setup decisions). When they're empty, the integration loop runs but contributes zero — exactly what we observe.

Workaround

For your specific case (uniform-viscosity Stokes with traction on the internal interface), you have a couple of options:

• Re-mesh as two stacked boxes with a shared boundary, apply traction as an exterior natural BC on each side, and rely on the geometric coincidence to act as the interface. Awkward but sound.
• Wait for the mesh-side fix (see below).

There isn't currently a tested UW3 workflow for traction on an interior interface — that's the territory of the "Interior penalty / Nitsche for internal boundaries" planning item, which acknowledges Nitsche-style enforcement is hard for internal facets. Natural BCs (just adding a known traction) should in principle be cleaner than Nitsche enforcement, but they're blocked here by the same label-population issue plus PETSc's interior-facet integration limitations.

Next steps

The right fix is in `BoxInternalBoundary`'s mesh-setup code (`meshing/cartesian.py`):

1. Expand each boundary label to include its closure (vertices bordering the labeled facets), matching what `StructuredQuadBox` does.
2. Populate `All_Boundaries` and `Null_Boundary` consolidated labels via `markBoundaryFaces` or equivalent.

Once the exterior labels are populated correctly, `add_natural_bc("Top")` etc. should work the same as on a vanilla mesh. Whether the internal facet case also works at that point is a separate question — PETSc's standard `DM_BC_NATURAL_FIELD` machinery handles facets with a single adjacent cell, and internal facets have two. That may need a custom assembly path (`BdIntegral`-style) for traction on internal interfaces.

For now, this needs to be picked up as a mesh-side fix; I'll let the team decide whether it's a small follow-up here or its own issue. Either way, no blocker on PETSc — the diagnosis is purely on the UW3 mesh-construction side.

Underworld development team with AI support from Claude Code

[lmoresi]

@NengLu — proposed fix in PR #163, ready for your review.

What the PR does

Fixes the label-population layer of this issue. `_dm_unstack_bcs` (used by `BoxInternalBoundary` to set up boundary labels after the gmsh import) was silently destroying `Null_Boundary` and `All_Boundaries` and not expanding the named labels' closure. After the fix, label populations on `BoxInternalBoundary` match `StructuredQuadBox` structurally:

label	SQBox	BIB pre-fix	BIB post-fix
`Top`	63	32	65
`All_Boundaries`	128	0	128
`Null_Boundary`	1089	0	1089

The same fix benefits the `mesh_adapt_meshVar` adaptation path.

Important caveat

This PR is necessary but may not be sufficient for your full case. After applying it, your reproducer still shows zero velocity in my local testing, suggesting the solver-side natural BC integration on `BoxInternalBoundary` has a deeper issue — likely related to recent normal-scaling changes (`Gamma_N` → `Gamma_P1`). The label fix is a strict prerequisite for any of that to work, so it's worth landing on its own.

Could you verify?

If you can pull `bugfix/bib-boundary-labels-162` and re-run your reproducer (the one in the issue body), the diagnostic we'd love to see is:

• Does the SNES initial F-norm now show a non-trivial value (not `< 1.e-11`)? That would confirm the natural BC contribution is at least being computed.
• If you set the natural BC on an exterior boundary instead of `Internal` (e.g. `Top`), does that work?

Either result tells us where the remaining work is. Once you confirm the label fix doesn't break your existing scripts — and we have a clearer picture of the residual issue — we'll merge the PR and close this issue, then track any follow-up under a new (more narrowly scoped) issue.

Note: GitHub couldn't add you as a formal reviewer (you're not a collaborator on the repo), so I'm @-mentioning here instead. Comment-back works the same way for our purposes.

Underworld development team with AI support from Claude Code

[NengLu]

Hi Louis,
Thank you. However, the test is still failing on the bugfix/bib-boundary-labels-162 branch:

1. The SNES initial F-norm shows an incorrect value (< 1.e-11, and u ≈ 2.5e-19 everywhere)
2. Setting the natural boundary condition on an exterior boundary (either Top or Bottom) causes the solver not to converge, resulting in u approaching zero everywhere.

for wall in ("Left", "Right", "Bottom"):
    stokes.add_essential_bc((0.0, 0.0), wall)
stokes.add_natural_bc(sympy.Matrix([1.0, 0.0]), "Top")

solver output:

  0 SNES Function norm < 1.e-11
  Nonlinear Solver_61_ solve converged due to CONVERGED_FNORM_ABS iterations 0
  0 SNES Function norm < 1.e-11
  Nonlinear Solver_61_ solve did not converge due to DIVERGED_LINE_SEARCH iterations 0
SNES solve diverged after 0 iterations: DIVERGED_LINE_SEARCH The solution vector may not have been updated. Use <solver>.petsc_options.setValue('snes_converged_reason', None) and <solver>.petsc_options.setValue('snes_monitor', None) to investigate.

[lmoresi]

@NengLu — update on PR #163: there was a second layer to the bug, now also fixed in the same PR.

What changed

The label-population fix was correct but not enough. There was a second issue: the Mesh constructor builds a consolidated `UW_Boundaries` label (the one the solver actually registers natural BCs against) at construction time — before `BoxInternalBoundary` runs its post-construction `_dm_unstack_bcs` call. So `UW_Boundaries` had only the auxiliary entries and none of the named boundaries (Top, Bottom, Left, Right, Internal). Every `add_natural_bc(...)` on a BIB mesh was being silently registered against an empty consolidated label.

Both layers are now fixed inside `_dm_unstack_bcs` (rebuild `UW_Boundaries` from the named labels at the end of the function).

Your reproducer

	Initial F-norm	u just above/below interface
Pre-fix	< 1e-11 (no BC contribution)	2.5e-19 (zero)
Post-fix	0.131 (real BC contribution)	0.0565 (matches your FEniCS reference)

The solver now converges in 1 SNES iteration on your exact case, and the magnitude matches your FEniCS analytical 5.65310750e-02 to roundoff. Sign is opposite (UW3 reports -0.0565); that's likely an internal-boundary normal-orientation convention rather than anything physically wrong, but worth your eye.

Could you verify?

If you can pull the updated `bugfix/bib-boundary-labels-162` and re-run your script:

1. Does the magnitude match the FEniCS reference? (I expect yes.)
2. Do any of your existing scripts regress? The fix touches the consolidated label rebuild, and Louis flagged that recent normal-scaling changes (`Gamma_N` → `Gamma_P1`) may have left some old scripts depending on subtly different label state.
3. The sign discrepancy (your FEniCS gives +0.0565, UW3 now gives -0.0565). If your FEniCS sign is unambiguously the physically expected one, we may want to track the orientation question as a separate issue — the magnitude part of [BUG] The natural boundary condition is not working on the internal boundary #162 is fixed regardless.

Once you confirm, we'll merge and close #162.

Underworld development team with AI support from Claude Code

[NengLu]

1. Yes!
2. So far, the old scripts show no new regression issues.
3. This is related to an orientation definition problem, but it's not a major issue.

Good news: The new fix also kind of makes the FSSA_internalBC work for ALEIB (only tested for using Gamma here, using Gamma_P1 or projected Normal term may need further testing).

It is okay for the merge! Thanks!

[lmoresi]

Fixed in #163 (merged to development). Two-layer bug: _dm_unstack_bcs was destroying both the auxiliary Null_Boundary/All_Boundaries labels and (separately) leaving the consolidated UW_Boundaries label without entries for the named boundaries. Solver registers natural BCs against UW_Boundaries, so on BoxInternalBoundary every add_natural_bc(...) was a silent no-op. Fix verified by @NengLu against his FEniCS reference: |u_x| = 0.0565 matches analytical reference to roundoff. Bonus: also unblocked FSSA_internalBC for ALEIB. Closing manually since the PR merged to development not main.