How to mark the boundary of a subdomain

Dear all,

I would like to mark the boundary of the subdomain, which is not the boundary of the entire domain. How can I achieve that?

I have tried below but it does not work

parameters['reorder_dofs_serial'] = False
x0, y0 = 10, 10 ### computational domain size

center_1 = np.array([-3.5, -4]) ### center of the hole domain
radius_1 = 0.5 ### radius of the hole domain

center_2 = np.array([3.5, 4])
radius_2 = 0.5


domain = Rectangle(Point(-x0,-y0),Point(x0,y0))
hole_1 = Circle(Point(center_1[0],center_1[1]),radius_1)
hole_2 = Circle(Point(center_2[0],center_2[1]),radius_2)

domain.set_subdomain(1, hole_1)
domain.set_subdomain(2, hole_2)

resol = 150 ### resolution of the mesh
mesh = generate_mesh(domain, resol) ### the FEM mesh
subdomains = MeshFunction('size_t',mesh,mesh.topology().dim(),mesh.domains()) ### subdomains data
n = FacetNormal(mesh) ### outward normal unit vector

V = FunctionSpace(mesh, "P", 1) ### Langrage base function is used


dx=Measure('dx',subdomain_data=subdomains)
submesh_hole_1 = SubMesh(mesh, subdomains,1)
submesh_hole_2 = SubMesh(mesh, subdomains,2)


############################################################################
### Define two subdomain boundaries separately
############################################################################
class interior_boundaries_1(SubDomain):
    def inside(self,x,on_boundary):
        return abs(sqrt((x[0]-center_1[0])**2+(x[1]-center_1[1])**2))-radius_1 <= mesh.hmin()

class interior_boundaries_2(SubDomain):
    def inside(self,x,on_boundary):
        return abs(sqrt((x[0]-center_2[0])**2+(x[1]-center_2[1])**2))-radius_2 <= mesh.hmin()



boundaries=MeshFunction('size_t',mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
interior_boundaries_1().mark(boundaries, 1)
interior_boundaries_2().mark(boundaries, 2)
ds=Measure('ds',subdomain_data=boundaries)

print(assemble(project(Constant(1), V) * ds(1))) ### expect pi but got 0.0

Actually, I was thinking to get the global index of the edges of the subdomain boundaries, then marking the boundaries manually, however, it does not recognize “parent_facet_indices”, while “parent_cell_indices” and “parent_vertex_indices” work well.

Any help is much appreciated!

Best,
Alice

If you save your boundaries to file and visualize them with Paraview, you will see that you have marked all interior facets of the mesh within the given radius.

File("mf.pvd") << boundaries

To me it is unclear if you only want to mark the interface between the interior circle and the rest of the domain, or every facet inside the domain. Please be more specific as to what you would like to achieve.

Additionally, as the facets are interior to the mesh, you need to use a dS measure, i.e

dS = Measure('dS', domain=mesh, subdomain_data=boundaries)

Dear dokken,

Thank you for your kind reply. I want to mark the boundary of the two holes, then do the integral over the hole boundary. I think the bracket position was wrong in the code… sorry about that.

Best,
Alice

See for instance: Plotting subdomains/internal boundaries - #2 by dokken