Edges in 3D not initialised in parallel
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.
Comments (9)
-
reporter -
reporter 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 callDistributedMeshTools::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()
inmesh.init()
. -
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.
-
reporter 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. -
Another reason to not silently perform a parallel operation is that local init becomes, unnecessarily, a collective operation.
-
reporter What's a suggested fix? I think the status quo is bad.
-
Throw an error if the requested entities haven't been computed.
-
reporter I've got a branch https://bitbucket.org/fenics-project/dolfin/branch/chris/fix-global-index which should fix this issue.
- 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. - Set global_size(dim) to zero unless set.
- Adds a method
Mesh::init_global(dim)
to create global indices on demand - otherwise you have to create aFunctionSpace
or something.
I'll create a pull request.
- Returns
-
reporter - changed status to resolved
- Log in to comment
I've added a test, in a branch
chris/add-3d-edge-mpi-test
which picks this up.