Hi, I have the following code, which solves Poisson problem and prints solution in a cell and it neighbours also building cells adjacency matrix which is needed later on. The solving part maps to firedrake just fine, but I'm unsure about what to do with loops over cells and their neighbours. Here the code: u = TrialFunction(V) v = TestFunction(V) f = Constant(1) a = inner(nabla_grad(u), nabla_grad(v))*dx L = f*v*dx # Compute solution u = Function(V) solve(a == L, u, bcs) tdim = mesh.topology().dim() mesh.init(tdim - 1, tdim) cell_neighbors = {cell.index(): sum((filter(lambda ci: ci != cell.index(), facet.entities(tdim)) for facet in facets(cell)), []) for cell in cells(mesh)} print cell_neighbors print mesh.cells() for cell in cells(mesh): point = cell.midpoint() #print cell.index(), list(cells(mesh))[cell_neighbors[cell.index()][0]] nhs = [u(Cell(mesh, cx).midpoint()) for cx in cell_neighbors[cell.index()]] print u(point), nhs Sincerely, George