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()

[kinnala]

Seems that you have an invisible node in the middle:

This plot can be produced by mesh.draw(node_numbering=True).show().

[duhd1993]

I found a very old thread in gmsh mailing list. They claim it has been fixed. But apparently it's not. Maybe it's worthwhile to filter out unused points in from_meshio?
http://geuz.org/pipermail/gmsh/2007/002566.html

[kinnala]

Mesh.is_valid() reveals this. I think Mesh.remove_unused() or Mesh.remove_duplicates() would be good when not automatic.

[duhd1993]

@kinnala Thanks! Are you going to implement remove_unused/remove_duplicates? The solution in gmsh is to assign a physical group to the entity to keep. They will only work when writing a file, adding unnecessary read/write cost.

[kinnala]

Not today, but I added an issue #1063 to remind myself (or someone else) to do it. I've been lately following a release cycle of 6 months unless there is something drastic.

[duhd1993]

Thanks!