I’m trying to get the neighbors of some cells in a dG mesh.
Neighbor in this case refers to cells sharing a face with the original cell.

In legacy fenics, you could just do something like

geodim = mesh.topology().dim()
mesh.init(geodim-1, geodim)
# subdomain is a dolfin.cpp.mesh.MeshFunctionSizet function indicating with 1/0 whether a cell does/does not belong to a subdomain of the mesh
for cell in df.cells(mesh):
if subdomain[cell.index()] == 1:
for face in df.facets(cell):
facet_cell_indices = face.entities(geodim)
for facet_cell_index in facet_cell_indices:
# here I can work with the neighboring cells since I have their cell index.

In fenicsx I haven’t found a way to iterate over cells of a mesh. Currently I’m working directly with dofs, but I cannot find a connection between the dofs of different cells. Is there a way to iterate over cells like in legacy fenics or get (at least one) dof of a neighboring cell if I have a dof inside a cell?

Any help would be much appreciated, thanks in advance!

I’m now struggling to get the mesh tags (or rather a certain face) from just a given number of dofs instead of a geometrical property.

In your example you got some faces by the geometrical property x[0] close to 0.5.

If I have only given some dofs I can identify the cells and the facets to that face, but how can I get the corresponding neighboring cells to a a facet without the ft meshtag map?/How can I create the ft map just being given the dofs of the nods?

example:

mesh = dolfinx.mesh.create_unit_interval(MPI.COMM_WORLD, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))

i.e. I have the unit interval with 10 cells (11 faces) and 2 dofs per cell. If I now only have the array

dofs_subdomain = np.array([2,3])

I can get from that, that the corresponding facet ids are 1 and 2, but how do I get the cell id from the other cells with facets 1 and 2?

How are you determining what dof indices you want to access? Are you assuming that the dofs are ordered left to right? Usually you cannot make that assumption (as meshes are unstructured and can be ran in both serial and parallel).

It was my understanding that mesh.topology.connectivity(mesh.topology.dim-1,mesh.topology.dim) gives me an adjacency list with of faces as nodes and the respective cells as values. mesh.topology.connectivity(mesh.topology.dim-1,mesh.topology.dim).offsets then gives me the indizes where the cell ids to a specific facet id in the mesh.topology.connectivity(mesh.topology.dim-1,mesh.topology.dim).array start and end.

Then you would need to invert the function-space dofmap (i.e. V.dofmap.cell_dofs(i) gives you want dofs in the function space is connected to cell i. You would have to loop over this, make the inverse map for each cell, and then map that to the facets of the cell.

What do you want to use the neighbouring cells for?

I don’t quite get how the snippet above connects to the dofs of your problem. In general there is no 1-1 mapping between a cell, facet or vertex and a degree of freedom of a function space. Please provide a minimal reproducible example.

import dolfinx
import numpy as np
from mpi4py import MPI
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 4)
V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))
mesh.topology.create_entities(mesh.topology.dim -1)
f_map = mesh.topology.index_map(mesh.topology.dim -1)
num_facets = f_map.size_local + f_map.num_ghosts
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
c_map = mesh.topology.index_map(mesh.topology.dim)
num_cells = c_map.size_local + c_map.num_ghosts
mesh.topology.create_connectivity(mesh.topology.dim, mesh.topology.dim-1)
c_to_f = mesh.topology.connectivity(mesh.topology.dim, mesh.topology.dim-1)
# Marking the dofs of two cells
dofs_submesh = np.zeros(V.dofmap.list.array.size, dtype=np.int32)
dofs_submesh[V.dofmap.cell_dofs(1)] = 1
dofs_submesh[V.dofmap.cell_dofs(17)] = 1
# New array for old cells + neighbors
neighbors = dofs_submesh.copy()
# This is the snippet from above
def get_cell_id(facet_id):
conn = mesh.topology.connectivity(mesh.topology.dim-1,mesh.topology.dim) # offsets, # array
offset_start = conn.offsets[facet_id]
offset_end = conn.offsets[facet_id + 1]
cell_ids = conn.array[offset_start:offset_end]
return cell_ids
# Collect all neighboring cells
cell_ids = set()
for cell in range(num_cells):
if 1 in dofs_subdomain[V.dofmap.cell_dofs(cell)]:
for facet in c_to_f.links(cell):
facet_cell_ids = get_cell_id(facet)
for facet_cell_id in facet_cell_ids:
cell_ids.add(facet_cell_id)
# Mark all new cells
for cell in cell_ids:
neighbors[V.dofmap.cell_dofs(cell)] = 1

If I then create a function with value 1 in the marked cells and 0 otherwise, the function is 1 in the two original cells and it’s face neighbors and 0 elsewhere: