apparent bug in assembling with mixed space for test functions

Question

Here is a simple piece of code that seems to indicate a bug. I want to compute a bilinear form acting on a test function u from Lagrange P1 on the unit interval, and a pair of test functions (v, c) from a mixed function space consisting of piecewise constants for v and a real number for c. The bilinear form has two parts:

  a0 = u.dx(0) * v * dx  # integral of u'*v over the interval
  a1 = u * c * Expression("1.-x[0]") * ds  # value of u*c at the point x=0

If the mesh has n intervals, then the first bilinear form should assemble into the first n rows of the (n+1) x (n+1) stiffness matrix, and the second into the final row. The problem is that in fact the second bilinear form assembles into the wrong row. Therefore if the full matrix is assembled from a0 + a1 it has a row of all zeros at the end, so is of course singular.

I see this both with version 1.3.0 and a development version installed on March 26.

from dolfin import *
mesh = UnitIntervalMesh(4)
V1 = FunctionSpace(mesh, 'CG', 1)
V2 = FunctionSpace(mesh, 'DG', 0)
R = FunctionSpace(mesh, 'R', 0)
u = TrialFunction(V1)
X = MixedFunctionSpace([V2, R])
(v, c) = TestFunctions(X)
# this bilinear form does not involve the real test function c,
# so one row should be identically zero
a0 = u.dx(0) * v * dx
print assemble(a0).array()
# this bilinear form only involves the real test function c
# so the corresponding row should have a nonzero
a1 = u * c * Expression("1 - x[0]") * ds
print assemble(a1).array()

Finally, here is the output.

[[ 0.  0.  0.  1. -1.]
 [ 0.  1. -1.  0.  0.]
 [ 0.  0.  1. -1.  0.]
 [ 1. -1.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]]
[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  1.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]]

Note that the zero row of the first array is the 5th row, but the nonzero row of the
2nd array is the third, not the fifth row.

Answer 1

This bug is fixed in https://bitbucket.org/fenics-project/dolfin/commits/3c0d0be1b571, and is now in the master branch of the DOLFIN development version. The fix will be in the DOLFIN 1.5 release.

Answer 2

Hi, I added a few diagnostics calls after you code

print 'V2'
print V2.dofmap().tabulate_all_coordinates(mesh)

print 'R'
print R.dofmap().tabulate_all_coordinates(mesh)

print 'X'
print X.dofmap().tabulate_all_coordinates(mesh)
print X.sub(0).dofmap().dofs()
print X.sub(1).dofmap().dofs()

These indicate, that somehow 0.375 dof of V2 is lost and nobody owns the last dof.
This explains the last row of zeros in the first matrix. Also the dof of R is 3rd, hence the nonzero third row of second matrix.

V2
[ 0.125 0.375 0.625 0.875]
R
[ 0.875]
X
[ 0.125 0.625 0.875 0.875 0. ]
[0 1 2 3]
[2]

This seems to be a bug in DofMapBuilder::reorder_local as (on dolfin 1.3.0) setting

parameters['reorder_dofs_serial'] = False 

fixes the problem.

[[-1.  1.  0.  0.  0.]
 [ 0. -1.  1.  0.  0.]
 [ 0.  0. -1.  1.  0.]
 [ 0.  0.  0. -1.  1.]
 [ 0.  0.  0.  0.  0.]]
[[ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 0.  0.  0.  0.  0.]
 [ 1.  0.  0.  0.  0.]]
V2
[ 0.125  0.375  0.625  0.875]
R
[ 0.875]
X
[ 0.125  0.375  0.625  0.875  0.875]
[0 1 2 3]
[4]