fix: enforce now handles rows with no nonzeros (fixes #1195)

[gaoflow]

Fix: enforce now handles rows with no nonzeros

Fixes #1195.

Problem

skfem.utils.enforce zeros the rows of A belonging to the Dirichlet
dofs D by constructing a flat array of positions into A.data. The
previous implementation:

idx = np.ones(count.sum(), dtype=np.int32)
idx[np.cumsum(count)[:-1]] -= count[:-1]
idx = np.repeat(start, count) + np.cumsum(idx) - 1

only works when every row in D has at least one stored nonzero. If any
count[i] == 0 — which happens whenever an enforced row has no in-bulk
neighbours, e.g. a MeshTri().refined(1) p-Laplacian Hessian whose
surface dofs have all-zero gradient contributions — then
np.cumsum(count)[:-1] contains duplicate indices, the corresponding
entries in idx are decremented twice, and the final idx runs past
the end of A.data and raises IndexError.

Repro from #1195 (a p-Laplacian Hessian on a 1-refined mesh, where the
surface dofs end up with empty rows):

import numpy as np, skfem as fem
from skfem.helpers import dot, grad

@fem.BilinearForm
def hessian(u, v, p):
    w = p['x']; gradw = grad(w); tmp = dot(gradw, gradw)
    plap = (4-2)*tmp**((4-4)*0.5) * dot(grad(u), gradw)*dot(gradw, grad(v))
    lap  = tmp**((4-2)*0.5)*dot(grad(u), grad(v))
    return lap + plap

mesh = fem.MeshTri().refined(1)
basis = fem.Basis(mesh, fem.ElementTriP1())
x = basis.project(lambda x: 0.02 * np.pi**2 * np.sin(np.pi*x[0]) * np.sin(np.pi*x[1]))
x[mesh.boundary_nodes()] = 0.
A = fem.asm(hessian, basis, x=basis.interpolate(x), p=4, rho=1.0)
fem.utils.enforce(A, D=mesh.boundary_nodes())
# IndexError: index 27 is out of bounds for axis 0 with size 27

Different from the closed #1043 (where the entire matrix is zero and
the result wouldn't be invertible) — here only a subset of the enforced
rows are empty and the resulting enforced system is perfectly usable.

Fix

Rewrite the index construction to handle empty rows correctly:

offsets = (np.arange(count.sum())
           - np.repeat(np.cumsum(count) - count, count))
idx = np.repeat(start, count) + offsets

Both np.repeat(start, count) and the offset arithmetic contribute
nothing when count[i] == 0, so the result is well-defined regardless
of how many zero-rows are in D. Also short-circuit when
count.sum() == 0 (e.g. D is empty) to avoid scipy warnings about
zero-length fancy indexing.

Verification

Adds tests/test_utils.py::test_enforce_handles_empty_rows. The test
sets up the failing scenario explicitly (a BilinearForm whose
coefficient is identically zero, which makes every boundary row empty)
and asserts:

• assert (counts == 0).any() — the test really triggers the broken
  code-path,
• enforced rows end up with diag == 1.0 and no off-diagonal nonzeros,
• non-enforced rows are bitwise-equal to the input,
• with a right-hand side, bout[D] == x[D] and other entries unchanged.

Local CI (matches .github/workflows/main.yml):

$ flake8 skfem
$ sphinx-build -a -b doctest docs docs/_build
... 188 tests, 0 failures in tests ...
$ pytest tests/test_utils.py tests/test_basis.py
======================= 87 passed, 104 warnings in 5.19s =======================

Scope

Only skfem/utils.py (production fix) and tests/test_utils.py
(regression test). No public API change — enforce's signature and
return contract are unchanged; the fix is purely the internal flat-index
construction.

[kinnala]

Thank you for the contribution!