system_assembler(): Bilinear and linear forms do not have same exterior facet subdomains in SystemAssembler.

[Nico Schlömer]

I'm getting warnings of the form

*** Warning: Bilinear and linear forms do not have same exterior facet subdomains in SystemAssembler. Taking subdomains from bilinear form

which I had overlooked in the deluge that is the application output. At other points in FEniCS, such inconsistencies are treated as errors (e.g., inconsistencies in dx). Shouldn't the same be the case here? (At least I would have profited from it.)

Here is an MWE that highlights how separate assembling of lhs and rhs works correctly, assemble_system() does not:

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)
ds_subdomain = Measure('ds')[submesh_boundaries]
#ds_subdomain = Measure('ds', domain=submesh)[submesh_boundaries]
#ds_subdomain = Measure('ds', domain=submesh, domain_data=submesh_boundaries)

u = TrialFunction(Vsub)
v = TestFunction(Vsub)
n = FacetNormal(submesh)
alpha = 100.0
beta = 100.0
F = dot(grad(u), grad(v)) * dx \
    - Constant(alpha) * v * ds_subdomain(0) \
    - Constant(beta) * v * ds_subdomain(1) \
    - f * v * dx

bc = DirichletBC(Vsub, 0.0, left)

# works:
#A = assemble(lhs(F))
#b = assemble(rhs(F))
# yields warning and produces the wrong result:
A, b = assemble_system(lhs(F), rhs(F))

bc.apply(A, b)

sol = Function(Vsub)

solver = KrylovSolver('cg', 'amg')
solver.parameters['relative_tolerance'] = 1.0e-12
solver.parameters['absolute_tolerance'] = 0.0
solver.solve(A, sol.vector(), b)
plot(sol)
interactive()

[Jan Blechta]

Duplicate of #78.