DirichletBC has not a proper way checking that tensors are compatible

Description: More careful compatibility checks in DirichletBC will be introduced by 670c1f385c27d5ce64e9123114baa33f945761a4 when solving ~~#717~~. But the problem is that in combination with ~~#719~~ user codes are prone to subtle bugs:

# Create tensors on one place
def create_data(A=None):
    "This function creates data used in the tests below"
    mesh = UnitSquareMesh(4, 4)
    V = FunctionSpace(mesh, "Lagrange", 1)
    u = TrialFunction(V)
    v = TestFunction(V)
    if A is None:
        A = assemble(u*v*dx)
    else:
        assemble(u*v*dx, tensor=A)
    b = assemble(v*dx)
    x = Vector()
    return A, x, b

# Create bc on other place from different dofmap!!!
def test_p31_box_2():
    mesh = UnitSquareMesh(4, 4)
    V = FunctionSpace(mesh, "CG", 1)
    A, x, b = create_data()
    u_0 = Expression("sin(x[0])", degree=1)
    bc = DirichletBC(V, u_0, "x[0] > 0.5 && on_boundary")
    u = Function(V)

    bc.apply(A, b)
    bc.apply(u.vector())

Although incorrect, this will worl fine in sequential and will produce random strange errors in parallel due to ~~#719~~. In fact, this is wrong code from test/unit/python/book/test_chapter10.py which was causing random failures for a long time and was not understood until ~~#717~~ was being resolved.

Possible resolution: One possibility would be that tensors would carry std::size_t GenericTensor::_index_map_id which would be compared to function_space->dofmap()->index_map()->id() in DirichletBC.

Comments (7)