Implementation of `GenericVector.__getitem__()` is wrong in parallel

[Jan Blechta]

from dolfin import *
mesh = UnitSquareMesh(3,3)
V = FunctionSpace(mesh, 'CG', 1)
u = Function(V)
x = u.vector()

pn = MPI.process_number()
n0, n1 = x.local_range()

# All local dofs including ghosts
dofs = V.dofmap().dof_to_vertex_map(mesh) + n0
print 'Process %i has %i dofs including ghosts.'%(pn, dofs.size)

# Set values [0,1,2,3,...]
for i in range(x.size()):
    x[i] = float(i)

# Only some values are printed because global vector
# of dofs.size is created - this is wrong.
print x[dofs].array()
x.update_ghost_values()
print x[dofs].array()

# This deadlocks because of global vector
# constructor waits for other processes.
if pn==0:
    print x[dofs].array()

[Jan Blechta]

Both demonstarted issues are solved by pull request 32.

[Johan Hake]

The index interface of GenericVector in Python is poorly tested in parallel. That is because it was written to work for serial cases, and we have not decided how it should work work in parallel. This issue illustrate this, and your pull request suggests one behavior:

Local vectors are returned when a accessing a sub-vector through the get{item,slice} interface in parallel.

I think that makes sense. Are there use cases where one would like it to return a global vector?

[Martin Sandve Alnæs]

Local makes sense, but it must be consistent with the rest of the interface. A typical use case would be to get a particular set of dof indices and use them for indexing. For example dofs with certain relations to topology and geometry and subdomains.

[Jan Blechta]

One could appreciate to get a global sub-vector corresponding to subspace.

Well, I didn't think it through a much. But if it passes tests we could perhapse consider it as a bugfix of current state which is obviously supposed to return local vector but has an implementation error.

[Johan Hake]

The GenericVector.array interface has always been local only. However, it has not been obvious that slices should return a local vector.

That said, I think it is a sane decision. For now a user need to make sure that the dofs resides on the local process.

Using a dofmap from a subspace a user can get the correct indices for what ever sub space he wants. The same can be acquired for dofs related to topology and geometry and subdomains. Creating nifty help functions for this in the interface is another issue.

[Jan Blechta]

Already fixed. Was not marked as resolved automatically from some reason :/