How to mark the boundary of a subdomain

Legacy DOLFIN
1

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

2

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

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

4

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

5

Hi there,
I saw some posts related to marking boundary conditions. Right now I try to mark the outer surface of a cylinder, but when I open it in Paraview, it does not show that a subdomain is actually allocated to this surface. Probably it is related to what is already mentioned over here, but I don’t seem to fix it.

from fenics import *
import mshr
from mshr import Cylinder

# Custom geometry definition
L = 2
R = 1
tol= 1e-2

mesh_cyl = Cylinder(Point(0, 0, 0), Point(0, 0, L), R, R)
mesh = mshr.generate_mesh(mesh_cyl, 50)

# Define subdomains
class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], 2)

class Surface(SubDomain):
    def inside(self, x, on_boundary):
        r = sqrt((x[0] - 1)**2 + (x[1] - 1)**2)
        return on_boundary and near(r, tol)

# Create MeshFunction for subdomains
subdomains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
subdomains.set_all(0)

# Mark subdomains
top = Top()
surface = Surface()

top.mark(subdomains, 1)
surface.mark(subdomains, 2)

# Write mesh and subdomains to XDMF file
xdmf_file = XDMFFile("mesh.xdmf")
xdmf_file.write(mesh)
xdmf_file.write(subdomains)
xdmf_file.close()

Or should it be better to make use of another code you have provided? If yes, what changes should be made beneath to the 2D case?

from dolfin import *
import numpy as np

def phi(x):
    return (x[0] - 0.5)**2 + (x[1]-0.5)**2 - 0.2**2

class Omega(SubDomain):
    def inside(self, x, on_boundary):
        return phi(x) > 0


mesh = UnitSquareMesh(30, 30)
domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)

omega = Omega() 
omega.mark(domains, 1)  
dx_int = Measure("dx", domain=mesh, subdomain_data=domains)

print(assemble(Constant(1)*dx_int(0)), assemble(Constant(1)*dx_int(1)))
File("domains.pvd") << domains


# Interior facets between the two domains
int_boundary = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
tdim = mesh.topology().dim()
mesh.init(tdim-1, tdim)
facet_to_cell = mesh.topology()(tdim-1, tdim)
num_facets = facet_to_cell.size()
domain_values = domains.array()
facet_values = int_boundary.array()
for facet in range(num_facets):
    # Check if interior
    cells = facet_to_cell(facet) 
    if len(cells == 2):
        # Check if facet is on the boundary between two cells with different markers
        values = np.unique(domain_values[cells])
        if len(values) == 2:
            facet_values[facet] = 1        
    else:
        continue

dS_int = Measure("dS", domain=mesh, subdomain_data=int_boundary) 

print(assemble(Constant(1)*dS_int(0)), assemble(Constant(1)*dS_int(1)))
File("facets.pvd") << int_boundary

I would like to hear from you.
Kind regards