Edges in 3D not initialised in parallel

[Chris Richardson]

There seems to be a bug in parallel in 3D for Edges.

from dolfin import *

mesh = UnitCubeMesh(5,5,5)
for dim in range(4):
    mesh.init(dim)
    print dim, mesh.size(dim), mesh.size_global(dim)

reports:

0 128 216
0 125 216
1 612 612
2 869 1650
3 384 750
1 591 591
2 833 1650
3 366 750

on two processes. Clearly the size(1) and size_global(1) should not be the same.

[Chris Richardson]

I've added a test, in a branch chris/add-3d-edge-mpi-test which picks this up.

[Chris Richardson]

OK, the problem seems to be that, although you can call mesh.init(1) to initialise the Edges, it does not do any global numbering, unless you also call DistributedMeshTools::number_entities(mesh, 1) (which is unavailable in Python).

As it happens, Edges in 3D are a special case, because vertices, cells and facets are always initialised globally in parallel.

My suggestion is that we should always call DistributedMeshTools::number_entities() in mesh.init().

[Prof Garth Wells]

Not sure that's a good solution. Potentially expensive calls should be made explicit.

Maybe we need separate calls for local and global, and we can also wrap them all behind another function.

[Chris Richardson]

Normally I'd agree, but the current situation leads to very non-intuitive behaviour. In parallel you already have numbered vertices, facets and cells globally. The only overhead would be when calling mesh.init(1) in parallel - but reasonably enough, as that is a fairly clear signal from the user. Also, it is hard (impossible?) from Python at the moment.

[Prof Garth Wells]

Another reason to not silently perform a parallel operation is that local init becomes, unnecessarily, a collective operation.

[Chris Richardson]

What's a suggested fix? I think the status quo is bad.

[Prof Garth Wells]

Throw an error if the requested entities haven't been computed.

[Chris Richardson]

I've got a branch https://bitbucket.org/fenics-project/dolfin/branch/chris/fix-global-index which should fix this issue.

1. Returns std::numeric_limits<std::size_t>::max() for unset global index (instead of crashing). This was in the comment, but not implemented until now.
2. Set global_size(dim) to zero unless set.
3. Adds a method Mesh::init_global(dim) to create global indices on demand - otherwise you have to create a FunctionSpace or something.

I'll create a pull request.

[Chris Richardson]

Fixed by https://bitbucket.org/fenics-project/dolfin/commits/91a037123d4f6453117951859a218adc48721050