Hi all,

I believe the simplest way to find the cells which share a facet is:

mesh.interior_facets.facet_cell

Cheers,

David

On Sat, 30 Apr 2016 at 19:09 Homolya, Miklós <m.homolya14@imperial.ac.uk> wrote:

Hi,

To evaluate a field at the midpoint of all cells, you can just interpolate onto a DG0 field. For example:

midpoints = interpolate(u, FunctionSpace(mesh, 'DG', 0))

Details: http://firedrakeproject.org/interpolation.html#ufl-expressions

Probably the simplest way to find cells that share a facet is to create a RT1 (RTCF1 on non-simplices) function space on the mesh, and inspect its cell_node_list attribute. It will be a #cells x #facets/cell matrix, if any two row contains the same number, that means they share a facet. The entries in midpoints.dat.data are ordered the same way as the rows of cell_node_list.

Regards,
Miklós

________________________________________
From: firedrake-bounces@imperial.ac.uk <firedrake-bounces@imperial.ac.uk> on behalf of George Ovchinnikov <lives9@gmail.com>
Sent: 30 April 2016 15:54:10
To: firedrake
Subject: [firedrake] Converting FENICS code to Firedrake (iterations over       mesh cells)

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

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake

_______________________________________________
firedrake mailing list
firedrake@imperial.ac.uk
https://mailman.ic.ac.uk/mailman/listinfo/firedrake