Change values of a Function only where a CellFunction has been marked.

Question

So I'm trying to simply set a function equal to zero where I have set a CellFunction equal to zero:

from fenics import *
from numpy  import where

mesh      = UnitCubeMesh(10,10,10)
Q         = FunctionSpace(mesh, 'CG', 1)
f         = Function(Q)
c         = CellFunction('size_t', mesh, 0)
c[-1]     = 1
dofmap    = Q.dofmap()
cells     = where(c.array() == 1)[0]
dofs      = []
for i in cells:
    dofs.extend(dofmap.cell_dofs(i))
dofs = list(set(dofs))

print MPI.rank(mpi_comm_world()), "dofs:", dofs

f.vector()[dofs] = 1.0

This works in serial, but in parallel I get a

IndexError: expected indices to be in [0..333]
IndexError: expected indices to be in [0..333]
IndexError: expected indices to be in [0..350]

I'm trying to find a local-to-global mapping, which is what I think I need to relate the dofmap.cell_dofs back to local vertices, right?

Comments

update - this works for with dolfin 1.4.0, but not with dolfin 1.5.0.

Answer 1

Hi, there has been changes to dof indexing in 1.5.0 version, see here. The following works for me with 1.5.0 and 1.6.0dev

from fenics import *

mesh = UnitSquareMesh(20, 20)
V = FunctionSpace(mesh, 'CG', 2)

cell_f = CellFunction('size_t', mesh, 0)
AutoSubDomain(lambda x: x[0] > x[1]).mark(cell_f, 1)
domain_cells = SubsetIterator(cell_f, 1)

dofmap = V.dofmap()
# Dofs in domain
dofs = sum((dofmap.cell_dofs(cell.index()).tolist()
           for cell in domain_cells), [])
# Get unique dofs in local numbering that the process sees. Some might not be owned
dofs = set(dofs)
# Where in global vector are dofs owned by the process
my_first, my_last = dofmap.ownership_range()
# Keep only owned ones, ie those whose global index is in ownership range
dofs = filter(lambda dof: my_first <= dofmap.local_to_global_index(dof) < my_last,
              dofs)
# Assign
f = Function(V)
f.vector()[dofs] = 1

plot(f)
interactive()