Hi,
I had the same problem. Here is an example that might help you:
from dolfin import *
from mshr import *
from math import pi
# -------------------- #
# Create mesh #
# -------------------- #
rectangle = Rectangle(Point(-1., -1.), Point(1., 1.))
R = 0.5
origin = Point(0.,0.)
circle = Circle(origin, R, segments=32)
domain = rectangle
domain.set_subdomain(1, circle)
domain.set_subdomain(2, rectangle - circle)
mesh = generate_mesh(domain, 15) # 15 is the resolution
# Create subdomains
subdomains = MeshFunction("size_t", mesh, 2, mesh.domains())
# -------------------- #
# Create boundaries #
# -------------------- #
class Left(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] + 1.) < DOLFIN_EPS
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[0] - 1.) < DOLFIN_EPS
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] + 1.) < DOLFIN_EPS
class Top(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and abs(x[1] - 1.) < DOLFIN_EPS
boundaries = FacetFunction("size_t", mesh, 0 ) # this index 0 is an alternative to the command boundaries.set_all(0)
Bottom().mark(boundaries, 1)
Left().mark(boundaries, 2)
Right().mark(boundaries, 2)
Top().mark(boundaries, 3)
# ---------------------------- #
# Create boundary interface #
# ---------------------------- #
# Option 1:
for f in facets(mesh):
domains = []
for c in cells(f):
domains.append(subdomains[c])
domains = list(set(domains)) # remove duplicates
if len(domains) > 1: # it is on the interface
boundaries[f] = 4
# Option 2:
#for f in facets(mesh): # For all edges in the mesh
# p0 = Vertex(mesh, f.entities(0)[0]) # save the two ending points p0 and p1
# p1 = Vertex(mesh, f.entities(0)[1])
# if near(p0.x(0)*p0.x(0) + p0.x(1)*p0.x(1), R*R, 0.05) and near(p1.x(0)*p1.x(0) + p1.x(1)*p1.x(1), R*R, 0.05): # check if the points lie on the circle - if yes, put a label on this edge
# boundaries[f] = 4
# --------------------------------- #
# Chech that everything is correct #
# --------------------------------- #
# Define a function of ones on the mesh
V = FunctionSpace(mesh, "CG", 1)
f = Function(V)
f.vector()[:] = 1.
# dS on the inteface = 2*pi*R
dS = dS(subdomain_data=boundaries)
form4 = f*dS(4) # integration on the whole facet marked with number 4
print form4
ans = assemble(form4)
print ans, "vs", 2*pi*R
# dS on an external boundary = 0
form1 = f*dS(1)
ans = assemble(form1)
print ans, "vs", 2., "~>", "dS is non zero only in the boundaries' interface!"
# ds on an external boundary = 2
ds = ds(subdomain_data=boundaries)
form1 = f*ds(1) # integration on the external boundary marked with number 1
ans = assemble(form1)
print ans, "vs", 2.