assembly of mixed spaces gives wrong results

[Johannes Ring]

Here @dnarnold reported the following:

Here is a simple piece of code high-lighting a bug which results in a stiffness matrix being completely wrong. It computes 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 first reported this in FEniCS Q&A for version 1.3.0. It persists in version 1.4.0 and the latest development version.

 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.

As an answer to the posting in Q&A MiroK added some diagnostic output that indicates the problems is in DofMapBuilder::reorder_local, and found that setting

parameters['reorder_dofs_serial'] = False

works around the problem.

[Martin Sandve Alnæs]

This code now gives:

martinal@martinal-mc:~/dev/fenics/dolfin/test/unit/python$ python tmp.py 
Traceback (most recent call last):
  File "tmp.py", line 7, in <module>
    X = MixedFunctionSpace([V2, R])
  File "/home/martinal/opt/fenics/dorsal-dev-141022/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 575, in __init__
    constrained_domain=constrained_domain)
  File "/home/martinal/opt/fenics/dorsal-dev-141022/lib/python2.7/site-packages/dolfin/functions/functionspace.py", line 168, in __init__
    dolfin_dofmap  = cpp.DofMap(ufc_dofmap, mesh)
  File "/home/martinal/opt/fenics/dorsal-dev-141022/lib/python2.7/site-packages/dolfin/cpp/fem.py", line 655, in __init__
    _fem.DofMap_swiginit(self,_fem.new_DofMap(*args))
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     [e-mail]
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to complete call to function build().
*** Reason:  Assertion MPI::sum(mesh.mpi_comm(), (std::size_t) dofmap._local_ownership_size) == dofmap._global_dimension failed.
*** Where:   This error was encountered inside ../../dolfin/fem/DofMapBuilder.cpp (line 184).
*** Process: unknown
*** 
*** DOLFIN version: 1.4.0+
*** Git changeset:  46798da7e9fa26c10c0b2662eed09729a00f9b6c
*** -------------------------------------------------------------------------

[Martin Sandve Alnæs]

Setting parameters['reorder_dofs_serial'] = False still makes the code work.

Changing the Real subspace to DG0 and keeping reordering on gives this message:

  File "/home/martinal/opt/fenics/dorsal-dev-141022/lib/python2.7/site-packages/dolfin/fem/assembling.py", line 203, in assemble
    assembler.assemble(tensor, dolfin_form)
RuntimeError: *** Error: All dofmaps must have same block size (for now)

[Martin Sandve Alnæs]

@garth-wells is this something you know about and plan to fix (considering the "for now" comment in the error)?

[Martin Sandve Alnæs]

Related to https://bitbucket.org/fenics-project/dolfin/issue/352/assembly-of-rectangular-matrices-broken

[Martin Sandve Alnæs]

Assigning to Garth because of the relation between this and https://bitbucket.org/fenics-project/dolfin/issue/352/assembly-of-rectangular-matrices-broken which Garth has said to fix.

[Prof Garth Wells]

Looking into this. Seems to be an issue with DG + Real elements. Below code throws an error (at least an error is better than a wrong result).

from dolfin import *
mesh = UnitIntervalMesh(4)
V2 = FunctionSpace(mesh, 'DG', 0)
R = FunctionSpace(mesh, 'R', 0)
X = MixedFunctionSpace([V2, R])

[Jan Blechta]

@garth-wells, with your failing example I observe using debugger that block size is 2. So my suspicion is that DofMapBuilder::compute_blocksize does not work as expected in this case.

[Prof Garth Wells]

@blechta , yes looks like the block size is not computed correctly. If I force bs=1, the resulting matrices appear correct. I'll figure out what's wrong with the block size.

[Prof Garth Wells]

Check for global dofs before computing dofmap block size.

Fixes #331. The 'actual' bug is in FFC/UFC, see https://bitbucket.org/fenics-project/ffc/issue/61/.

→ <<cset 3c0d0be1b571>>

[Prof Garth Wells]

Removing milestone: 1.5 (automated comment)