Mesh topology connections return non-existing cell numbers

[Tormod Landet]

This could very possibly be my fault for not initializing the mesh correctly, but the following code works in serial and reports non-existing cell numbers when run on two processors. I expected either the correct results or some error message

from dolfin import *

parameters['ghost_mode'] = 'shared_vertex'

comm = MPI.comm_world
rank = MPI.rank(comm)

mesh = UnitSquareMesh(2, 2)
mesh.init(1, 2)

# Print cell neighbours for each facet and count the number
for facet in facets(mesh):
    fidx = facet.index()
    cnbs = mesh.topology()(1, 2)(fidx)
    print('Rank %d facet %d has connected cells %s' % (rank, fidx, cnbs))

In parallel I get results like

Rank 0 facet 0 has connected cells [37763008]
Rank 0 facet 1 has connected cells [0 3105283032]
Rank 0 facet 2 has connected cells [32520]
Rank 0 facet 3 has connected cells [0 4]
Rank 0 facet 4 has connected cells [1 2]
...

[Tormod Landet]

With latest master version I sometimes get reasonable results, but on most runs I get random results for at least some facets. I have tried all types of mesh.init() combinations just to be sure.

[Francesco Ballarin]

I think this only occurs with pybind11 wrappings, or at least I run your snippet a bunch of times with swig wrappings and everything was working correctly.

[Tormod Landet]

You are probably right, I just assumed for whatever reason that this was related to the ghosting, but I get the same results without ghosting the mesh (no errors in serial, though). I am indeed using the pybind11 wrappers

[Jan Blechta]

    py::class_<dolfin::MeshConnectivity, std::shared_ptr<dolfin::MeshConnectivity>>
      (m, "MeshConnectivity", "DOLFIN MeshConnectivity object")
      .def("__call__", [](const dolfin::MeshConnectivity& self, std::size_t i)
           {
             return Eigen::Map<const Eigen::Matrix<unsigned int, Eigen::Dynamic, 1>>(self(i), self.size(i));})

is referencing plain ptr

std::size_t MeshConnectivity::size(std::size_t entity) const
const unsigned int* MeshConnectivity::operator() (std::size_t entity) const

Add py::return_value_policy::reference_internal.

[Tormod Landet]

Will do. I'm writing a test for it now first, fix and test will follow soon

[Jan Blechta]

You can try testing reference counts: https://bitbucket.org/fenics-project/dolfin/pull-requests/429/fix-pybind11-lifetime-bug-in-meshtopology/diff#Ltest/unit/python/mesh/test_mesh.pyT639.

[Tormod Landet]

Fixed by d3a9b1b