# How to integrate on the surface of a subdomain (without making the surface its own subdomain)

Hello All,

I am working with a multi-domain mesh where the subdomains are defined by a MeshFunction that is marked by XDMF cell attributes. Because of this, I cannot easily describe the domain surfaces as their own subdomains.

I want to integrate over the surface of a subdomain using a ‘dS’ or ‘ds’ measure I have tried the following, using a simple MeshFunction to make the subdomains:

``````from dolfin import *

mesh = UnitSquareMesh(20,20)

#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# toy example creating a MeshFunction sub_domains, similar to what I read in
# as a mesh attribute from an XDMF mesh
tol = 1e-14
class Right(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0] > 0.5 + tol
sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim())
sub_domains.set_all(0)
Right().mark(sub_domains,1)
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

x = SpatialCoordinate(mesh)
f = Constant(1)*x[1]

R_full = f*ds
print(assemble(R_full))
# output is as expected: 2.0

ds_custom = Measure('ds',domain=mesh,subdomain_data=sub_domains,subdomain_id=1)
dS_custom = Measure('dS',domain=mesh,subdomain_data=sub_domains,subdomain_id=1)

R_subdomain_s =f*ds_custom
print(assemble(R_subdomain_s))
# output is 0.0, should be 1

R_subdomain_S =f*dS_custom
print(assemble(R_subdomain_S))
# output is 0.03249999999999996, should be 1
``````

What is wrong with how I am defining these measures? Is there another way FEniCS has to integrate on the surfaces of subdomains?

You issue is in the definition of the `MeshFunction` as you define it as a function over cells.

You should change it to

``````sub_domains = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
sub_domains.set_all(0)
Right().mark(sub_domains,1)
File("mf.pvd") << sub_domains
``````

and you can visualize the markers on facets in for instance Paraview.
Using this marker will get you something close to 1 for

However, I disagree with the statement that

should be 1, as `dS` is the integration measure over internal facets (those that are not on the boundary).
As no internal markers has the value 1, the integral will be 0

Thank you for your reply! I now see the how the distinction between ‘ds’ and ‘dS’ will work.

Unfortunately, in the application I am working on, I am constructing the MeshFunction with mesh attributes from my mesh XDMF file, with the domain markers living on the cells. So when I define the MeshFunction as you suggest, the cell IDs jumble up with the facet IDs and I get a garbage output.

I see now I will have to change how I am dealing with the domain data to associate it with facets instead. Thanks again for your assistance!