Boundary condition lifting does not respect domain makers

[jorgensd]

Lifting for problems with subdomains have unclear behavior, as it only uses the kernel for the first integral:
https://github.com/FEniCS/dolfinx/blob/master/cpp/dolfinx/fem/assemble_vector_impl.h#L174

Consider a 3x1 unit square mesh with quadrilateral, where the cells are marked 1,0,2, from left to right.
Assume that we have a boundary condition on the left and right hand side of the unit square.
Assembling each of the domains 1 and 2 separately, and applying the boundary condition can only reproduce the expected behavior if there is only one cell integral present. See the minimal example below:

import numpy as np
import ufl
from mpi4py import MPI
from petsc4py import PETSc
import dolfinx
import dolfinx.io
def left_side(x):
    return x[0] < 1/3
def right_side(x):
    return x[0] > 2/3
mesh = dolfinx.UnitSquareMesh(
    MPI.COMM_WORLD, 3, 1, dolfinx.cpp.mesh.CellType.quadrilateral)
V = dolfinx.FunctionSpace(mesh, ("Lagrange", 1))

tdim = mesh.topology.dim
num_cells = mesh.topology.index_map(tdim).size_local
cell_midpoints = dolfinx.cpp.mesh.midpoints(mesh, tdim,
                                            range(num_cells))
values = []
for cell_index in range(num_cells):
    if left_side(cell_midpoints[cell_index]):
        values.append(1)
    elif right_side(cell_midpoints[cell_index]):
        values.append(2)
    else:
        values.append(0)
ct = dolfinx.mesh.MeshTags(mesh, mesh.topology.dim,
                           range(num_cells),
                           np.array(values, dtype=np.intc))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
x = ufl.SpatialCoordinate(mesh)
dx = ufl.Measure("dx", domain=mesh, subdomain_data=ct)
a = x[0]*ufl.inner(u, v)*dx(1) + 2*ufl.inner(u, v) * dx(2)
a1 = x[0]*ufl.inner(u, v)*dx(1)
a2 = 2*ufl.inner(u, v)*dx(2)
l = 2*ufl.inner(x[0], v)*dx(1) +\
    ufl.inner(x[1], v)*dx(2)
l1 = 2*ufl.inner(x[0], v)*dx(1)
l2 = ufl.inner(x[1], v)*dx(2)


def right(x):
    return np.isclose(x[0], 1)


def left(x):
    return np.isclose(x[0], 0)


left_facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, left)
right_facets = dolfinx.mesh.locate_entities_boundary(
    mesh, mesh.topology.dim-1, right)
boundary_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, np.hstack([left_facets, right_facets]))
u1 = dolfinx.Function(V)
with u1.vector.localForm() as u_local:
    u_local.set(1)
bc12 = dolfinx.DirichletBC(u1, boundary_dofs)

left_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, left_facets)
bc1 = dolfinx.DirichletBC(u1, left_dofs)
right_dofs = dolfinx.fem.locate_dofs_topological(
    V, mesh.topology.dim-1, right_facets)

bc2 = dolfinx.DirichletBC(u1, right_dofs)


# Assemble vector for each case
L = dolfinx.fem.assemble_vector(l)
L1 = dolfinx.fem.assemble_vector(l1)
L2 = dolfinx.fem.assemble_vector(l2)
L12 = dolfinx.fem.assemble_vector(l)

dolfinx.fem.apply_lifting(L, [a], [[bc12]])
L.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
              mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(L, [bc12])

dolfinx.fem.apply_lifting(L12, [a1, a2], [[bc1], [bc2]])
L12.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(L12, [bc1, bc2])

dolfinx.fem.apply_lifting(L1, [a1], [[bc1]])
L1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
               mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(L1, [bc1])

dolfinx.fem.apply_lifting(L2, [a2], [[bc2]])
L2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
               mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(L2, [bc2])

Lbc1 = dolfinx.fem.assemble_vector(l)
dolfinx.fem.apply_lifting(Lbc1, [a], [[bc1]])
Lbc1.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                 mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(Lbc1, [bc1])

Lbc2 = dolfinx.fem.assemble_vector(l)
dolfinx.fem.apply_lifting(Lbc2, [a], [[bc2]])
Lbc2.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES,
                 mode=PETSc.ScatterMode.REVERSE)
dolfinx.fem.set_bc(Lbc2, [bc2])

print("-"*25, "L", "-"*25, "\n", L.array)
print("-"*25, "L1", "-"*25, "\n", L1.array)
print("-"*25, "L2", "-"*25, "\n", L2.array)
print("-"*25, "L1+L2", "-"*25, "\n", L1.array+L2.array)
print("-"*25, "L12", "-"*25, "\n", L12.array)
assert(np.allclose(L12.array, L1.array+L2.array))

Here, we would expect that the L array, which used the bc for both domains and the matrix for both domains, should reproduce the sum of each of the integrals, L1 and L2.
This is not the case:

------------------------- L ------------------------- 
 [1.         0.00462963 0.03240741 1.         0.03240741 0.03240741
 1.         1.        ]
------------------------- L1 ------------------------- 
 [0.         0.         0.         0.         0.03240741 0.03240741
 1.         1.        ]
------------------------- L2 ------------------------- 
 [ 1.00000000e+00 -2.77777778e-02  1.38777878e-17  1.00000000e+00
  0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00]
------------------------- L1+L2 ------------------------- 
 [ 1.00000000e+00 -2.77777778e-02  1.38777878e-17  1.00000000e+00
  3.24074074e-02  3.24074074e-02  1.00000000e+00  1.00000000e+00]
------------------------- L12 ------------------------- 
 [ 1.00000000e+00 -2.77777778e-02  1.38777878e-17  1.00000000e+00
  3.24074074e-02  3.24074074e-02  1.00000000e+00  1.00000000e+00]

[michalhabera]

Fixed via #1177