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.