Mshr ignores subdomains

I have previously used “mshr” to generate simple polygonal subdomains, and used it to generate a mesh that includes these subdomains. A few years later, the exact same code seems to now ignore the subdomains when generating the mesh. Are anyone having similar issues?

The following code “should” make a unit circle mesh with a subdomain perfectly aligned with a pentagon:

"""
N-sided polygon of "radius" R and center 'C' = [cx,cy]
"""
def polygon_array(N,C,R):
    n = np.array(range(N))
    pz = R*np.exp(2.0*n*1j*pi/N)
    P_array = []
    for i in n:
        P_array.append(Point( np.real(pz[i])+C[0], np.imag(pz[i])+C[1] ))
    return P_array

"""
Construct circle mesh with center 'center' and radius 'R' and fineness 'mesh_fineness'.
'inc_list' is a collection of boundaries for inclusions.
"""
def inclusion_mesh(center,R,mesh_fineness, bdry_fineness, inc_list):
    C = Polygon(polygon_array(bdry_fineness,center,R))
    k = 0
    for P in inc_list:
        k = k+1
        C.set_subdomain(k,P)
    mesh = generate_mesh(C,mesh_fineness)
    return mesh

fineness = 50
bdry_pts = 500
C = [.4,-.4]
R_pol = 3.0/14.0
N = 5
P = Polygon(polygon_array(N,C,R_pol))
mesh = inclusion_mesh([0.0,0.0],1.0,fineness,bdry_pts,[P]) 

When I use

subdomains = MeshFunction(“size_t”, mesh, 2, mesh.domains())
plot(subdomains)

It shows something that resembles a pentagon as the subdomain, but clearly the mesh generation has not taken the exact geometry into account, as the edges are very jagged.

example

Mshr is unfortunately no longer maintained, and users are adviced to move to other (external) meshing software, such as GMSH or pygmsh.

OK that is unfortunate. I prefer not to use too much external software. However I still find it strange that previously working code suddenly isn’t working anymore; I was mainly wondering if other people noticed the same.

The script is working fine, the problem comes from the matplotlib output which is not so great regarding resolution (don’t ask me why). If you export your plot to a high-res .png or a .svg file you will see that it is fine.

plt.savefig("fig.png", dpi=600)

2 Likes

Thanks, you are absolutely correct. Great!

I must admit, I am a bit embarrassed as I have used far too long to try and figure this out before posting. But I am also surprised that the default plot with matplotlib is of such bad quality.