CSG Geometry set_subdomain() code and docs disagree -- which is correct?

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()

Note: I previously this as a question -- https://fenicsproject.org/qa/12862/geometry-set_subdomain-code-and-docs-disagree-which-correct?show=12862

Comments (2)