Add `transfer_meshtags_to_submesh`

[jorgensd]

Describe new/missing feature

With the new co-dim 0 submesh support, it would be beneficial to have a function that transfers meshtags to the submesh (similar to what happens for mesh refinement). Interface (here written in Python for illustration purposes) is described below. It should be written in C++ as it involves alot of looping.

Suggested user interface

def transfer_meshtags_to_submesh(mesh, entity_tag, submesh, sub_vertex_to_parent, sub_cell_to_parent):
    """
    Transfer a meshtag from a parent mesh to a sub-mesh.
    """

    tdim = mesh.topology.dim
    cell_imap = mesh.topology.index_map(tdim)
    num_cells = cell_imap.size_local + cell_imap.num_ghosts
    mesh_to_submesh = np.full(num_cells, -1)
    mesh_to_submesh[sub_cell_to_parent] = np.arange(len(sub_cell_to_parent), dtype=np.int32)
    sub_vertex_to_parent = np.asarray(sub_vertex_to_parent)

    submesh.topology.create_connectivity(entity_tag.dim, 0)

    num_child_entities = submesh.topology.index_map(
        entity_tag.dim).size_local + submesh.topology.index_map(entity_tag.dim).num_ghosts
    submesh.topology.create_connectivity(submesh.topology.dim, entity_tag.dim)

    c_c_to_e = submesh.topology.connectivity(submesh.topology.dim, entity_tag.dim)
    c_e_to_v = submesh.topology.connectivity(entity_tag.dim, 0)

    child_markers = np.full(num_child_entities, 0, dtype=np.int32)

    mesh.topology.create_connectivity(entity_tag.dim, 0)
    mesh.topology.create_connectivity(entity_tag.dim, mesh.topology.dim)
    p_f_to_v = mesh.topology.connectivity(entity_tag.dim, 0)
    p_f_to_c = mesh.topology.connectivity(entity_tag.dim, mesh.topology.dim)
    for facet, value in zip(entity_tag.indices, entity_tag.values):
        facet_found = False
        for cell in p_f_to_c.links(facet):
            if facet_found:
                break
            if (child_cell := mesh_to_submesh[cell]) != -1:
                for child_facet in c_c_to_e.links(child_cell):
                    child_vertices = c_e_to_v.links(child_facet)
                    child_vertices_as_parent = sub_vertex_to_parent[child_vertices]
                    is_facet = np.isin(child_vertices_as_parent, p_f_to_v.links(facet)).all()
                    if is_facet:
                        child_markers[child_facet] = value
                        facet_found = True
    tags =  dolfinx.mesh.meshtags(submesh, entity_tag.dim,
                                 np.arange(num_child_entities, dtype=np.int32), child_markers)
    tags.name = entity_tag.name
    return tags

Minor edits: Some extra connectivity calls and variable renaming have taken place after testing this with @ottarph

[jorgensd]

A C++ implementation of this has been prototyped at: https://github.com/scientificcomputing/scifem/pull/73 Any feedback on how it should be altered to be part of DOLFINx is appreciated.

[dajuno]

Is this feature on the horizon?

[jorgensd]

Is this feature on the horizon?

As I haven't gotten any feedback, I haven't moved it to DOLFINx yet. You can use the scifem version, as scifem is on conda, spack and pip.