Vertex assembly does not work with mixed Lagrange x Real function space

[Ivan Yashchuk]

Assembly over vertices fails if one of the subspaces has dofs other than on vertices. Example:

from fenics import *

mesh = UnitSquareMesh(2, 2)
P1 = FiniteElement('P', mesh.ufl_cell(), 1)
R = FiniteElement('Real', mesh.ufl_cell(), 0)
W = FunctionSpace(mesh, MixedElement([P1, R]))

u, c = TrialFunctions(W)
v, d = TestFunctions(W)

Then assemble(u*dP) gives error. If lines 512-523 in dolfin/fem/Assembler.cpp are ignored, then it kinda works.

In [..]: assemble(u*dP).get_local()
Out[..]: array([ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.])

But the actual problem arrives if I try assemble(u*d*dP) or assemble(v*c*dP), the output is zeros in this case. I expect it to work the same way as ds and dx works.

In [..]: assemble(u*d*dP).array()
Out[..]:
array([[ 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.,  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.],
       [ 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.,  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.]])

The expected matrix:

array([[ 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.,  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.],
       [ 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.,  0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  0.]])

This is needed to have vertex constraints in 2D and 3D. Any tips on fixing that?

[Michal Habera]

Please note, that this type of "vertex" constraints do not always have correct mathematical background. In 2D (even more in 3D) there is no meaning for point evaluation of H^1 Sobolev functions. You need better regularity of the solution space. I wouldn't wonder if this casts some troubles also at the FEniCS level.

[Jan Blechta]

This is not a bug report but a feature request. Needs changes (or complete rewrite) of Assembler::assemble_vertices. The method only considers dofs associated with vertices:

632       // Get local dofs of the local vertex
633       dofmaps[i]->tabulate_entity_dofs(local_to_local_dofs[i], 0, local_vertex);

[Jan Blechta]

P1 space as well as R space is $C^0$ so evaluation at point has a good meaning.