scikit-fem: Singular matrix when mesh contains points not used by cells.

Expected behavior is that these points would be ignored (as is the case in other packages, e.g. SfePy, dolfinx). Observed behavior is assembly of a singular matrix and failure to solve.

Minimum code to reproduce is below, with one line to add/not add the offending point. You can inspect mesh.points and mesh.cells_dict['lines'] and mesh.cells_dict['triangles'] to see this point is in the points array but not used by any line or trinagle.

import gmsh
import meshio
import skfem
from skfem.helpers import dot, grad
from skfem.io.meshio import from_meshio

lc0 = .25
p_center = [.1234, .5678]
gmsh.initialize()
try:
    gmsh.model.add('mesh')
    # Define unit square area
    # Corners
    bottom_left = gmsh.model.geo.addPoint(0, 0, 0, lc0)
    bottom_right = gmsh.model.geo.addPoint(1, 0, 0, lc0)
    top_right = gmsh.model.geo.addPoint(1, 1, 0, lc0)
    top_left = gmsh.model.geo.addPoint(0, 1, 0, lc0)
    # Edges
    bottom = gmsh.model.geo.addLine(bottom_left, bottom_right)
    right = gmsh.model.geo.addLine(bottom_right, top_right)
    top = gmsh.model.geo.addLine(top_right, top_left)
    left = gmsh.model.geo.addLine(top_left, bottom_left)
    # Surface
    perimeter = gmsh.model.geo.addCurveLoop([bottom, right, top, left])
    surface = gmsh.model.geo.addPlaneSurface([perimeter])
#
#
# comment out the following line to remove the offending point
# and clear the warnings/failure 
#
    gmsh.model.geo.add_point(p_center[0], p_center[1], 0, lc0)
#
#
#
    gmsh.model.geo.synchronize()
    gmsh.model.mesh.generate(dim=2)
    gmsh.write('demo.msh')
finally:
    gmsh.finalize()

mesh = meshio.read('demo.msh')
fem_mesh = from_meshio(mesh)
basis = skfem.Basis(fem_mesh, skfem.ElementTriP1())

@skfem.BilinearForm
def laplace(u, v, _):
    return dot(grad(u), grad(v))

@skfem.LinearForm
def rhs(v, _):
    return 1.0 * v

A = laplace.assemble(basis)
b = rhs.assemble(basis)

D = fem_mesh.boundary_nodes()
A,b = skfem.enforce(A, b, D=D)  # warning generated here:
# site-packages\scipy\sparse\_index.py:125: SparseEfficiencyWarning: Changing the
# sparsity structure of a csr_matrix is expensive. lil_matrix is more efficient.
#    self._set_arrayXarray(i, j, x)
#

x = skfem.solve(A, b)  # another warning generated here:
# site-packages\scipy\sparse\linalg\dsolve\linsolve.py:206: MatrixRankWarning:
# Matrix is exactly singular
#    warn("Matrix is exactly singular", MatrixRankWarning)