The mshr documentation for CSGGeometry::set_subdomain() is very clear:
Define subdomain. This feature is 2D only. The subdomain is itself a
CSGGeometry and the corresponding cells in the resulting will be
marked with i If subdomains overlap, the latest added will take
precedence.
However, the implementation does the opposite. Example below shows that rectangle r1 is preserved even though r2 (added later) overlaps it. Is the bug in the docs or in the code?
from dolfin import *
from mshr import *
domain = Rectangle(Point(0,0), Point(10,10))
r1 = Rectangle(Point(1,1), Point(5,5))
r2 = Rectangle(Point(4,4), Point(6,6))
domain.set_subdomain(1, r1)
domain.set_subdomain(2, r2)
mesh = generate_mesh(domain, 1)
File("mesh.xml") << mesh
plot(mesh)
mf = MeshFunction('size_t', mesh, 2, mesh.domains())
File("mf.xml") << mf
plot(mf)
interactive()
