Ordering of cells via mesh.topology.original_cell_index

General
1

Hello,

I read in another post that mesh.topology.original_cell_index should return the prescribed ordering as imposed by the user before dolfinx changes order internally. In this little example, I do not find myself with the results, meaning that the ordering which I get in original_cell_idx is not always correct. The ordering which I need for further purposes is the one which scrolls among rows first, starting from the bottom one. Is there anything that I am missing here?

import dolfinx
import numpy as np
from mpi4py import MPI
import ufl

gdim = 2
shape = "quadrilateral"
degree = 1
N = 5

cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
X, Y = np.meshgrid(np.linspace(0, 1, N), np.linspace(0, 1, N))
X, Y = X.flatten(), Y.flatten()
coords = np.vstack([X, Y]).transpose()
items = [item for item in np.arange(N ** 2 - N - 1) if (item+1) % N != 0]
cells = [[i, i+1, i+N, i+N+1] for i in items]
cells = np.array(cells, dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, coords, domain)

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

original_cell_idx = mesh.topology.original_cell_index
2

I’m not really sure how you have tested the correctness of this mapping, as it gives the expected vertices for each cell for me:

import ufl
import numpy as np
import dolfinx
from mpi4py import MPI
N = 5

shape = "quadrilateral"
gdim = 2
degree = 1
cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))
X, Y = np.meshgrid(np.linspace(0, 1, N), np.linspace(0, 1, N))
X, Y = X.flatten(), Y.flatten()
coords = np.vstack([X, Y]).transpose()
items = [item for item in np.arange(N ** 2 - N - 1) if (item+1) % N != 0]
cells = [[i, i+1, i+N, i+N+1] for i in items]
cells = np.array(cells, dtype=np.int32)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, coords, domain)


num_cells_on_process = mesh.topology.index_map(mesh.topology.dim).size_local
geometry_indices = dolfinx.cpp.mesh.entities_to_geometry(mesh, mesh.topology.dim, np.arange(num_cells_on_process, dtype=np.int32), False)
for i, o_index in enumerate(mesh.topology.original_cell_index):
    print(f"Local cell {i}, input {o_index}", "Input vertices", coords[cells[o_index]], "Mesh vertices", mesh.geometry.x[geometry_indices[i]])
    assert(np.allclose(coords[cells[o_index]], mesh.geometry.x[geometry_indices[i]][:,:gdim]))

print(mesh.topology.original_cell_index)

gives:

Local cell 0, input 0 Input vertices [[0.   0.  ]
 [0.25 0.  ]
 [0.   0.25]
 [0.25 0.25]] Mesh vertices [[0.   0.   0.  ]
 [0.25 0.   0.  ]
 [0.   0.25 0.  ]
 [0.25 0.25 0.  ]]
Local cell 1, input 1 Input vertices [[0.25 0.  ]
 [0.5  0.  ]
 [0.25 0.25]
 [0.5  0.25]] Mesh vertices [[0.25 0.   0.  ]
 [0.5  0.   0.  ]
 [0.25 0.25 0.  ]
 [0.5  0.25 0.  ]]
Local cell 2, input 4 Input vertices [[0.   0.25]
 [0.25 0.25]
 [0.   0.5 ]
 [0.25 0.5 ]] Mesh vertices [[0.   0.25 0.  ]
 [0.25 0.25 0.  ]
 [0.   0.5  0.  ]
 [0.25 0.5  0.  ]]
Local cell 3, input 2 Input vertices [[0.5  0.  ]
 [0.75 0.  ]
 [0.5  0.25]
 [0.75 0.25]] Mesh vertices [[0.5  0.   0.  ]
 [0.75 0.   0.  ]
 [0.5  0.25 0.  ]
 [0.75 0.25 0.  ]]
Local cell 4, input 5 Input vertices [[0.25 0.25]
 [0.5  0.25]
 [0.25 0.5 ]
 [0.5  0.5 ]] Mesh vertices [[0.25 0.25 0.  ]
 [0.5  0.25 0.  ]
 [0.25 0.5  0.  ]
 [0.5  0.5  0.  ]]
Local cell 5, input 8 Input vertices [[0.   0.5 ]
 [0.25 0.5 ]
 [0.   0.75]
 [0.25 0.75]] Mesh vertices [[0.   0.5  0.  ]
 [0.25 0.5  0.  ]
 [0.   0.75 0.  ]
 [0.25 0.75 0.  ]]
Local cell 6, input 3 Input vertices [[0.75 0.  ]
 [1.   0.  ]
 [0.75 0.25]
 [1.   0.25]] Mesh vertices [[0.75 0.   0.  ]
 [1.   0.   0.  ]
 [0.75 0.25 0.  ]
 [1.   0.25 0.  ]]
Local cell 7, input 6 Input vertices [[0.5  0.25]
 [0.75 0.25]
 [0.5  0.5 ]
 [0.75 0.5 ]] Mesh vertices [[0.5  0.25 0.  ]
 [0.75 0.25 0.  ]
 [0.5  0.5  0.  ]
 [0.75 0.5  0.  ]]
Local cell 8, input 9 Input vertices [[0.25 0.5 ]
 [0.5  0.5 ]
 [0.25 0.75]
 [0.5  0.75]] Mesh vertices [[0.25 0.5  0.  ]
 [0.5  0.5  0.  ]
 [0.25 0.75 0.  ]
 [0.5  0.75 0.  ]]
Local cell 9, input 12 Input vertices [[0.   0.75]
 [0.25 0.75]
 [0.   1.  ]
 [0.25 1.  ]] Mesh vertices [[0.   0.75 0.  ]
 [0.25 0.75 0.  ]
 [0.   1.   0.  ]
 [0.25 1.   0.  ]]
Local cell 10, input 7 Input vertices [[0.75 0.25]
 [1.   0.25]
 [0.75 0.5 ]
 [1.   0.5 ]] Mesh vertices [[0.75 0.25 0.  ]
 [1.   0.25 0.  ]
 [0.75 0.5  0.  ]
 [1.   0.5  0.  ]]
Local cell 11, input 10 Input vertices [[0.5  0.5 ]
 [0.75 0.5 ]
 [0.5  0.75]
 [0.75 0.75]] Mesh vertices [[0.5  0.5  0.  ]
 [0.75 0.5  0.  ]
 [0.5  0.75 0.  ]
 [0.75 0.75 0.  ]]
Local cell 12, input 13 Input vertices [[0.25 0.75]
 [0.5  0.75]
 [0.25 1.  ]
 [0.5  1.  ]] Mesh vertices [[0.25 0.75 0.  ]
 [0.5  0.75 0.  ]
 [0.25 1.   0.  ]
 [0.5  1.   0.  ]]
Local cell 13, input 11 Input vertices [[0.75 0.5 ]
 [1.   0.5 ]
 [0.75 0.75]
 [1.   0.75]] Mesh vertices [[0.75 0.5  0.  ]
 [1.   0.5  0.  ]
 [0.75 0.75 0.  ]
 [1.   0.75 0.  ]]
Local cell 14, input 14 Input vertices [[0.5  0.75]
 [0.75 0.75]
 [0.5  1.  ]
 [0.75 1.  ]] Mesh vertices [[0.5  0.75 0.  ]
 [0.75 0.75 0.  ]
 [0.5  1.   0.  ]
 [0.75 1.   0.  ]]
Local cell 15, input 15 Input vertices [[0.75 0.75]
 [1.   0.75]
 [0.75 1.  ]
 [1.   1.  ]] Mesh vertices [[0.75 0.75 0.  ]
 [1.   0.75 0.  ]
 [0.75 1.   0.  ]
 [1.   1.   0.  ]]
[ 0  1  4  2  5  8  3  6  9 12  7 10 13 11 14 15]
3

Hi @dokken,

Thanks for your answer.

I was diving into it now and the problem is that the vertices (or at least the degrees of freedom as I see them numbered in Paraview, I know they may not coincide) are renumbered as well, hence the current code does not give the desired result. Is there any way in which I can fix this, to the best of your knowledge? Ideally, I need somehow to have an ordering of the cells.

import dolfinx
import numpy as np
from mpi4py import MPI
import ufl

def F(x, y, z, mod):
    theta = mod.pi
    return [mod.cos(theta) * x - mod.sin(theta) * y,
            mod.sin(theta) * x + mod.cos(theta) * y,
            z]

gdim = 2
shape = "quadrilateral"
degree = 1
N = 5

cell = ufl.Cell(shape, geometric_dimension=gdim)
domain = ufl.Mesh(ufl.VectorElement("Lagrange", cell, degree))

X, Y = np.meshgrid(np.linspace(0, 1, N), np.linspace(0, 1, N))
X, Y = X.flatten(), Y.flatten()
coords = np.stack([X, Y]).transpose()
items = [item for item in np.arange(N ** 2 - N - 1) if (item+1) % N != 0]
cells = [[i, i+1, i+N, i+N+1] for i in items]
cells = np.array(cells, dtype=np.int32)
print(cells)
mesh = dolfinx.mesh.create_mesh(MPI.COMM_WORLD, cells, coords, domain)

mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-2)
map_cells_to_vertices = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-2)
for cell in np.arange(len(cells)):
    print(f"Cell {mesh.topology.original_cell_index[cell]} has vertices {map_cells_to_vertices.links(cell)}")

Namely, I am expecting this code to give the same results when printing cells above and below, but this is not the case.