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
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
Dear George, On 30/04/16 15:54, George Ovchinnikov wrote:
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
David and Miklos have described how you can find this information in firedrake (although cell neighbours change depending on whether you mean neighbours through faces or vertices). However, can you explain what it is you intend to do with it? Perhaps there is a better way. Cheers, Lawrence
participants (3)
-
George Ovchinnikov
-
Homolya, Miklós
-
Lawrence Mitchell