Extracting Facet Dofs might lead to wrong results

Hello,

i think i found a bug that must be somehow related to the dofmap. I have a code, that extracts the dofs in a 3D Nedelec Function Space("N1curl") that belong to a facet for which facet.exterior() is true, hence it must lie on the boundary. If i set the regarding dofs to something ~0 with:

test= Function(V3D)
test.vector()[facet_dofs] = 1

the plot should show a vector that is parallel to the boundary (for Nedelec Elements of first kind), but however it does not. Here my full minimal code:

from dolfin import *


mesh = UnitCubeMesh(3,3,3)
mesh.init()
V3D = FunctionSpace(mesh,"N1curl", 1)

cell_number = 0

tdim = mesh.topology().dim()
dofmap = V3D.dofmap()
facet_dofs = [None]*len(dofmap.tabulate_facet_dofs(0))

cell = Cell(mesh,cell_number)
facets = cell.entities(tdim-1)
for local_facet_index, global_facet_index in enumerate(facets):
        if (Facet(mesh,global_facet_index).exterior()):
            local_facet_dofs = dofmap.tabulate_facet_dofs(local_facet_index)
            facet_dofs= dofmap.cell_dofs(cell_number)[local_facet_dofs]

facet_dofs = list(set(facet_dofs))

test= Function(V3D)
test.vector()[facet_dofs] = 1  

plot(test)
interactive()

It might be aswell an error that comes somehow from the interpolation of the plot (since he interpolates the solution to the nodes/vertices).

I guess the error is not related to Facet.exterior() since i can check if the facet is exterior "manually" by checking for the midpoint (and knowing where the boundary is) and get the same results.

The part with the facet_dofs from cell is taken from: http://fenicsproject.org/qa/2328/is-there-an-analog-of-dofmap-cell_dofs-cell-index-for-facets

Unfortunately i could not further reduce the problem nor figure out, where it really comes from. I apologize in advance if i am simply to stupid and this is not a bug.

Kind regards

Johann

Comments (5)