Singular when there is hole in model

[duhd1993]

linsolve.py:206: MatrixRankWarning: Matrix is exactly singular This is a minimal example using 2D elasticity. K is singular when I add a hole. I believe there is no mathematical difficulty to cause that. Even if I comment out holes parameters in add_polygon, I still get the error. I have to comment out the circle part as well. This makes me think it's because how from_meshio handles meshio objects. I'm still trying to debug. Any help is appreciated.

from skfem import *
from skfem.io import from_meshio
from skfem.helpers import *
import pygmsh
from skfem.visuals.matplotlib import draw, plot
from skfem.models.elasticity import linear_elasticity, linear_stress

def create_rect():
    with pygmsh.geo.Geometry() as geom:
        circle = geom.add_circle(
            x0=[0.5, 0.5, 0.0],
            radius=0.25,
            mesh_size=0.3,
            num_sections=4,
            make_surface=False,
        )
        rect = geom.add_polygon(
            [[0, 0, 0], [1, 0, 0], [1, 1, 0], [0, 1, 0]],
            mesh_size = 0.3,
            holes = [circle.curve_loop]
        )

        # geom.set_recombined_surfaces([rect.surface])
        mesh = geom.generate_mesh()
    return mesh

mesh = from_meshio(create_rect())

e = ElementTriP1()
ev = ElementVector(e)
ib = Basis(mesh, ev)

dofs = {
    'left' : ib.get_dofs(lambda x: x[0] == 0.0),
    'right': ib.get_dofs(lambda x: x[0] == 1.0),
}

E = 200e9
nu = 0.3
lam = E*nu/(1+nu)/(1-2*nu)
mu = E/2/(1+nu)
K = asm(linear_elasticity(lam, mu), ib)

u = ib.zeros()
u[dofs['right'].nodal['u^1']] = 0.3

I = ib.complement_dofs(dofs)

u = solve(*condense(K, x=u, I=I))
u1 = u[ib.get_dofs(elements=True).all(['u^1'])]

plot(mesh, u1, colorbar=True).show()