identify adjacent cells for set of cells

Question

Given a set C of cells of some mesh, I would like to find the layer of adjacent cells, i.e., all cells that share at least one vertex or edge with one of the cells in C. Apparently, iterating all cells in C and testing each with all cells of the mesh is too costly. Any suggestions how to achieve this efficiently?
Cheers, Martin

Answer 1

Not optimal this either, but doesn't check all cells in mesh. If C is quite small, this should do it.

from dolfin import *
from numpy import argwhere

mesh = UnitCubeMesh(4,4,4)

# Create meshfunction of cells
mf = CellFunction('size_t', mesh)
mf.set_all(0)
for i, c in enumerate(cells(mesh)):
    if c.midpoint()[0] < 0.5:
        mf[i] = 1

C = argwhere(mf.array()==1)

# Check adjacent cells for all cells in C
mesh.init(3,3)
for index in C:
    index = index[0]
    cell = Cell(mesh, index)
    for adj_cell in cell.entities(3):
        if mf[int(adj_cell)] == 0:
            mf[int(adj_cell)] = 2
plot(mf, interactive=True)

Comments

I guess this will do it for me. Thanks.

---

This variant is significantly faster.

tic()
mesh.init(3,3)
mf2 = CellFunction('size_t', mesh)
mf2.set_all(0)
for index in cells:
    index = index[0]
    cell = Cell(mesh, index)
    mf2.array()[cell.entities(3)] = 1
mf.array()[:] += mf2.array()
print toc()
plot(mf, interactive=True)

---

nice, thanks!

Answer 2

Try:

mesh = UnitSquareMesh(2,2)
for cell0 in cells(mesh):
    print cell0.index(), cell0.entities(0)
    for cell1 in cells(cell0):
        print "  ", cell1.index(), cell1.entities(0)

Comments

I just realized I did not answer your questions...

---

well, it's the intention that counts.

---

Besides, by that I learned that cells() can also be used with a single cell entity which comes in handy for determining the patch of a cell. Will employ it in another question ;)