MeshEntityIterator in DOLINx

dolfinx
1

Hi,

Is the MeshEntityIterator avaliable？For the new DOLFINx (v0.6.0), is there any functions that can be applied to do the entity iteration?

Thanks.

2

As DOLFINx is designed with functional programming in mind, you can easily make an iterator with the following:
entities = range(mesh.topology.index_map(dim).size_local) if you want to iterate over all entities of dimension dim owned by the current process.

If you could provide a minimal example of what entity information you want to obtain, it would be possible to Give a more specific advice

4

I find some functions in the source codes here https://github.com/FEniCS/dolfinx/blob/v0.6.0/cpp/dolfinx/fem/assemble_vector_impl.h

Screenshot from 2023-10-09 15-26-51
Screenshot from 2023-10-09 15-26-51532×693 45 KB

this function is what I want which enable me do do iteration over cells and access the DOFs in each cell. If I want to do the similar operation over edges, do you have some suggestions? Thanks

5

Just as @dokken wrote

You can iterate over entities of arbitrary dimension. If your problem has an edge topology of dimension dim_edge then the edge iterator would iterate over

entities = range(mesh.topology.index_map(dim_edge).size_local)

For example, a 3D problem would have dim = 3 and dim_edge = 1.

6

An example of how to get the dofs on the facets (or any sub entity of the facet) is shown below:

import dolfinx
import numpy as np
from mpi4py import MPI


mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10)

V = dolfinx.fem.VectorFunctionSpace(mesh, ("Lagrange", 2))
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: (x[0], 2*x[1]))

dim = mesh.topology.dim -1 # Use facets in example

# Compute cell->entity and entity->cell connectivity
mesh.topology.create_connectivity(dim, mesh.topology.dim)
mesh.topology.create_connectivity(mesh.topology.dim, dim)
dim_to_cell = mesh.topology.connectivity(dim, mesh.topology.dim)
cell_to_dim = mesh.topology.connectivity(mesh.topology.dim, dim)

# dof layout
dl = V.dofmap.dof_layout
assert (dl.block_size == V.dofmap.bs)
bs = dl.block_size

# Get dof layout for each cell
num_entities = mesh.topology.index_map(dim).size_local + mesh.topology.index_map(dim).num_ghosts
entity_dofs = []
for entity in range(num_entities):
    cells = dim_to_cell.links(entity) # Find all cells connected to entity
    entities = cell_to_dim.links(cells[0]) # Pick one cell to find facet local index in cell
    local_index = np.flatnonzero(entities==entity)
    closure_dofs = dl.entity_closure_dofs(dim, local_index) # Get dofs on facet for that cell
    dofs = V.dofmap.cell_dofs(cells[0])[closure_dofs]

    entity_dofs.append([]) # Store all dofs on facet (unrolled)
    for dof in dofs:
        for i in range(bs):
            entity_dofs[-1].append(dof*bs+i)


# Test by finding single facet

def one_facet(x):
    return (x[0]>=(0.9- 1e-10)) & np.isclose(x[1], 1)

facet = dolfinx.mesh.locate_entities(mesh, dim, one_facet )
assert len(facet) == 1
facet_dofs = entity_dofs[facet[0]]
x = V.tabulate_dof_coordinates()
print(u.x.array[facet_dofs], x[np.array(facet_dofs[0::2])//bs])
2 Likes