Getting neighbors/dofs of neighbor cells in dG mesh

[dergl]

Hello all,

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!

[dokken]

Consider the following minimal example:

import dolfinx
import numpy as np
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)
V = dolfinx.fem.FunctionSpace(mesh, ("DG", 1))
u = dolfinx.fem.Function(V)
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
facets = np.arange(num_facets, dtype=np.int32)
values = np.zeros_like(facets, dtype=np.int32)

def interface(x):
    return np.isclose(x[0], 0.5)

mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
values[dolfinx.mesh.locate_entities(mesh, mesh.topology.dim-1, interface)] = 1

ft = dolfinx.mesh.meshtags(mesh, mesh.topology.dim-1, facets, values)

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)
for cell in range(num_cells):
    for facet in c_to_f.links(cell):
        if ft.values[facet] == 1:
            print(cell, facet, V.dofmap.cell_dofs(cell))
            u.x.array[V.dofmap.cell_dofs(cell)] = 2

with dolfinx.io.VTXWriter(mesh.comm, "u.bp", [u], engine="BP4") as bp:
    bp.write(0.0)

giving:

[dergl]

Thanks, that helped a lot!

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?

[dokken]

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).

[dergl]

I have two methods:

• I take dofs of cells that are smaller than a certain threshold.
• Some geometrical property (as in your example)

I then want to take these cells and the neighboring cells (and at some point, also the neighbors of these neighbors)

[dergl]

Okay I think I solved it (if I understood mesh.topology.connectivity correctly, it appears to be working at least in 1d).

Assuming I’m given the facet_id and mesh as above, I can do the following to get the cell ids of the facet-neighbors:

conn = mesh.topology.connectivity(mesh.topology.dim-1,mesh.topology.dim)
offset_start = conn.offsets[facet_id]
offset_end = conn.offsets[facet_id + 1]
cell_ids = conn.array[offset_start:offset_end]

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.

[dokken]

@dergl

I take dofs of cells that are smaller than a certain threshold.

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?

[dergl]

I only need the dofs of the neighboring cells, nothing else.

The code in the comment I wrote above your reply seems to be working (in 1d and 2d at least), shouldn’t this work? Or did I miss something?

[dokken]

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.

[dergl]

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:

fct = dolfinx.fem.Function(V)
fct.x.array[np.where(neighbors==1)[0]] = 1

What I was interested in where the dofs/indices of the dofs, i.e.

np.where(neighbors==1)[0]

[dokken]

@dergl

if 1 in dofs_subdomain[V.dofmap.cell_dofs(cell)]:

dofs_subdomain is not defined in this script, so it is unclear to me how you are fetching those.

[dergl]

Sorry, it was a typo, I fixed it.