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.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 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?