Incorrect cell coordinate dofs on ghosted meshes

[Jan Blechta]

Following snippet run with -np 4

from dolfin import *
parameters["ghost_mode"] = "shared_vertex"
mesh = UnitSquareMesh(32, 32)
r = MPI.rank(mesh.mpi_comm())
if r == 0:
    c = Cell(mesh, 566)
    print(r, c.global_index(), c.get_coordinate_dofs())
if r == 3:
    c = Cell(mesh, 51)
    print(r, c.global_index(), c.get_coordinate_dofs())

gives output

0 1042 [ 0.3125   0.5      0.28125  0.5      0.3125   0.53125]
3 1042 [ 0.28125  0.5      0.3125   0.5      0.3125   0.53125]

and demonstrates that cell with global index 1042 gives wrong coordinate dofs. Coordinate dofs must be in the same order on all ranks as dictated by UFC. In fact value on rank 0 has permuted first and second vertex.

As a result tabulate_tensor functions are fed with rubbish and return rubbish.

[Jan Blechta]

Check if this fixes #835, #837.

[Jan Blechta]

@chris_richardson, this has something to do with ordering of mesh topology/geometry data. Any insight into the problem?

[Jan Blechta]

from dolfin import *
parameters["ghost_mode"] = "shared_vertex"
mesh = UnitSquareMesh(32, 32)
r = MPI.rank(mesh.mpi_comm())
if r == 0:
    c = Cell(mesh, 566)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), c.get_coordinate_dofs(), vertices_global)
    for v in vertices(c):
        print(r, v.global_index(), v.point().x(), v.point().y())
if r == 3:
    c = Cell(mesh, 51)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), c.get_coordinate_dofs(), vertices_global)
    for v in vertices(c):
        print(r, v.global_index(), v.point().x(), v.point().y())

shows that geometry data are correct. Cell::entities(0) needs fixing.

[Jan Blechta]

Partial fix in https://bitbucket.org/fenics-project/dolfin/branch/jan/fix-issue-840#diff. The problem is still that the mesh thinks it is ordered by the cached result while unordered ghost entities are computed later:

from dolfin import *
parameters["ghost_mode"] = "shared_vertex"
set_log_level(TRACE)
mesh = UnitSquareMesh(32, 32)
r = MPI.rank(mesh.mpi_comm())

# Bad ordering
if r == 0:
    c = Cell(mesh, 566)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), vertices_global)
if r == 3:
    c = Cell(mesh, 51)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), vertices_global)

assert not mesh.ordered()  # Fails!!!
mesh.order()
assert mesh.ordered()

# Good ordering
if r == 0:
    c = Cell(mesh, 566)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), vertices_global)
if r == 3:
    c = Cell(mesh, 51)
    vertices_global = [Vertex(mesh, v).global_index() for v in c.entities(0)]
    print(r, c.global_index(), vertices_global)

[Jan Blechta]

8c1c6a5cab20c526a2b0b75a3265cf868e9a59ad