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

1

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?

2

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

3

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!