skfem.enforce fails with zero matrix

[gatling-nrl]

mesh = skfem.MeshTri()
basis = skfem.CellBasis(mesh, skfem.ElementTriP1())

@skfem.BilinearForm
def form1(u, v, w):
    return w.e * dot(grad(v), grad(u))

A = form1.assemble(basis, e=np.zeros(basis.N))
# A[A==0] = 0   # bizarrely, uncomment this line to make the exception go away.
b = np.zeros(basis.N)
skfem.enforce(A, b, D=basis.get_dofs())

skfem\utils.py:375, in enforce(A, b, x, I, D, diag, overwrite)
    373 count = stop - start
    374 idx = np.ones(count.sum(), dtype=np.int64)
--> 375 idx[np.cumsum(count)[:-1]] -= count[:-1]
    376 idx = np.repeat(start, count) + np.cumsum(idx) - 1
    377 Aout.data[idx] = 0.

IndexError: index 0 is out of bounds for axis 0 with size 0

[gdmcbain]

I'm away from a terminal now but is it that mesh here is very small and all boundary so that all its degrees of freedom are enforced by the essential condition?

[gdmcbain]

No, I see that it is indeed only because A is a ‘zero matrix’, not because D.size == A.shape[0]; the behaviour is essentially unchanged by appending .refined() to the first line.

[gatling-nrl]

That's right. And i first noticed it on a much larger gmsh generated mesh. You might wonder how this came up, and it's because I'm extending some codes to support complex numbers, but when using those codes to run purely real things or purely imaginary things, you end up a zero matrix like this.

[gdmcbain]

Yes, I was wondering momentarily about applying essential conditions to degrees of freedom in the kernel of an operator, but I decided that that might not be up to enforce to question, the choice being left to the user.

A fix looks easy enough. The problem is that enforce is trying to zero rows that are already zero—or rather (and this explains why the hack A[A==0]=0 works) not there in the sparse representation of A.

[kinnala]

I don't think this should work. Such matrix cannot be invertible. Do you think we should add explicit zeros to the diagonal of all inputs? What is the use case?