Obtain high order mesh edge and face nodes?

mesh
1

Hi everyone,
I’m reading the dofmap paper Construction of Arbitrary Order Finite Element Degree-of-Freedom Maps on Polygonal and Polyhedral Cell Meshes | ACM Transactions on Mathematical Software and given a high order mesh, I’m wondering how I can obtain the edge and face nodes as in for example Fig 12 and Table 3? Here I can see edge #3 contain vertices 19, 18, 17. How can I get such info from the code?

I’m testing with a P2 quadrilateral mesh consisting of one P2 Lagrange element, and mesh.topology.connectivity and dolfinx.cpp.mesh.entities_to_geometry only give me the “linear” nodes. Any hints on which function I should use?

Thanks,
August

2

I made a function for this that was unfortunately not merged a while back: Extend and clean-up entities to geometry by jorgensd · Pull Request #1684 · FEniCS/dolfinx · GitHub

You should be able to do these computations in Python.
Maybe @mscroggs has a convenience script for this somewhere?

3

Thanks, that looks like some interesting code. I’ll have a look and see what I can decrypt!

4

Here’s a python variant of the function in the PR:

def entities_to_geometry(mesh, dim, entity_list):
    # See https://github.com/FEniCS/dolfinx/pull/1684

    layout = mesh.geometry.cmap.create_dof_layout()
    num_entity_dofs = layout.num_entity_closure_dofs(dim)
    xdofs = mesh.geometry.dofmap

    tdim = mesh.topology.dim
    mesh.topology.create_entities(dim)
    mesh.topology.create_connectivity(dim, tdim)
    mesh.topology.create_connectivity(tdim, dim)
    e_to_c = mesh.topology.connectivity(dim, tdim)
    c_to_e = mesh.topology.connectivity(tdim, dim)

    entity_geometry = np.empty((len(entity_list), num_entity_dofs))

    for i, idx in enumerate(entity_list):
        cell = e_to_c.links(idx)[0]
        cell_entities = c_to_e.links(cell)
        pos = np.nonzero(cell_entities == idx)[0]
        assert pos.size
        local_entity = cell_entities[pos]
        entity_dofs = layout.entity_closure_dofs(dim, local_entity)

        xc = xdofs.links(cell)
        entity_geometry[i] = xc[entity_dofs]

    return entity_geometry
2 Likes
5

This could probably be improved by using numba, and pure numpy arrays for e_to_c and c_to_e, xdofs etc. (if one pre-computes all the entity_closure dofs and store them in a 2D array).

However, as long as the code is sufficiently fast for you, no reason to change a winning recipe.