facet specification doesn't work on `ds(submesh)`

[Nico Schlömer]

When using the ds(submesh) notation, boundaries are not recognized as marked and Dolfin silently does the wrong thing.

MWE:

from dolfin import *

mesh = RectangleMesh(0.0, 0.0, 1.0, 1.0, 100, 100)


class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[1], (0.5, 0.7)) and between(x[0], (0.5, 1.0)))
obstacle = Obstacle()
subdomains = CellFunction('size_t', mesh)
subdomains.set_all(0)
obstacle.mark(subdomains, 1)

V = FunctionSpace(mesh, 'CG', 1)
f = project(Constant(1.0), V)

submesh = SubMesh(mesh, subdomains, 1)

Vsub = FunctionSpace(submesh, 'CG', 1)
v = TestFunction(Vsub)


class Left(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary \
            and x[0] < 0.5 + DOLFIN_EPS
left = Left()


class Right(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary \
            and x[0] > 1.0 - DOLFIN_EPS
right = Right()

submesh_boundaries = FacetFunction('size_t', submesh)
submesh_boundaries.set_all(0)
right.mark(submesh_boundaries, 1)

u = TrialFunction(Vsub)
v = TestFunction(Vsub)
n = FacetNormal(submesh)
alpha = 2.0
beta = 1.0
F = dot(grad(u), grad(v)) * dx(submesh) \
    - alpha * v * ds(submesh)(0) \
    - beta * v * ds(submesh)(1) \
    - Constant(0.0) * v * dx(submesh)

sol = Function(Vsub)
solve(lhs(F) == rhs(F), sol,
      bcs=DirichletBC(Vsub, Expression('x[0]*x[0]'), Left())
      )
plot(sol)
interactive()

[Martin Sandve Alnæs]

You're not using submesh_boundaries anywhere.

[Martin Sandve Alnæs]

User program bug.

[Nico Schlömer]

I thought I do via right. What would be the correct code?

[Nico Schlömer]

More importantly, still, I think that it is not good that these things fail silently. If the user makes a mistake, that fact should be made very clear.

[Martin Sandve Alnæs]

If you create an object some place in your program and don't use it anywhere, there's no way we can discover that mistake.

Here are the only three lines referencing submesh_boundaries:

submesh_boundaries = FacetFunction('size_t', submesh)
submesh_boundaries.set_all(0)
right.mark(submesh_boundaries, 1)

This amounts to creating an array, setting all values to 0, then setting some values to 1. Then you throw it away, by not using it anywhere.

[Nico Schlömer]

I'd been under the impression that marking subdomains does something under the hood that influences how ds operates.

I suppose the (a?) correct way of doing things would be to

submesh_boundaries = FacetFunction('size_t', submesh)
submesh_boundaries.set_all(0)
right.mark(submesh_boundaries, 1)
ds_subdomain = Measure('ds')[submesh_boundaries]

and then use ds_subdomain for defining the form. Is that correct?

[Martin Sandve Alnæs]

That's one way. Also these should work:

ds = Measure("ds", domain=submesh)[submesh_boundaries]
ds = Measure("ds", domain=submesh, domain_data=submesh_boundaries)
a = f*ds(1)

but we're pondering renaming domain_data to subdomain_data which would be more descriptive.

[Martin Sandve Alnæs]

"marking subdomains does something under the hood that influences how ds operates."

That has never been the case.

[Nico Schlömer]

Alright, thanks for the clarification. My form now looks like

F = dot(grad(u), grad(v)) * dx(submesh) \
    - Constant(alpha) * v * ds_subdomain(0) \
    - Constant(beta) * v * ds_subdomain(1) \
    - f * v * dx(submesh)

where I have to use dx(submesh) instead of dx_subdomain as in

subdomain = CellFunction('size_t', submesh)
subdomain.set_all(0)
dx_subdomain = Measure('dx')[subdomain]

because f lives on the supermesh (cf. the bugs you recently fixed).