# Facet map from submesh to parent mesh

**URL:** <https://fenicsproject.discourse.group/t/facet-map-from-submesh-to-parent-mesh/15118>\
**Category:** dolfinx\
**Created:** [July 1, 2024, 3:37pm UTC](https://fenicsproject.discourse.group/t/facet-map-from-submesh-to-parent-mesh/15118 "2024-07-01T15:37:57Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ordinary-slim](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ordinary-slim/32/4821_2.png) [@ordinary-slim](https://fenicsproject.discourse.group/u/ordinary-slim)\
**Post date:** [July 1, 2024, 3:37pm UTC](https://fenicsproject.discourse.group/t/facet-map-from-submesh-to-parent-mesh/15118/1 "2024-07-01T15:37:57Z")

</div>

I’m trying to build a mapping between facets in a submesh and a parent mesh. I tried emulating code from `create_submesh` which returns the same correspondences for cells and nodes. This Python code works in serial but doesn’t work in parallel:

```auto
from mpi4py import MPI
from dolfinx import mesh, fem, cpp, io
import numpy as np

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

def main():
    points_side = 16
    big_mesh = mesh.create_unit_square(MPI.COMM_WORLD,
                                        points_side,
                                        points_side,
                                        )
    cdim = big_mesh.topology.dim
    v_big = fem.functionspace(big_mesh, ("CG", 1),)
    dg0_big = fem.functionspace(big_mesh, ("Discontinuous Lagrange", 0),)
    left_els = fem.locate_dofs_geometrical(dg0_big, lambda x : x[0] <= 0.5 )
    submesh, subcell_map, subvertex_map, subgeom_map = mesh.create_submesh(big_mesh,cdim,left_els)
    left_els_facets = mesh.compute_incident_entities(big_mesh.topology,
                                                     subcell_map,
                                                     cdim,
                                                     cdim-1)
    facet_index_map, subfacet_map = cpp.common.create_sub_index_map(big_mesh.topology.index_map(cdim-1),
                                                    left_els_facets,
                                                    False)
    submesh.topology.set_index_map(cdim-1,facet_index_map)
    bfacets = mesh.locate_entities_boundary(submesh,1,lambda x : np.ones(x.shape[1]))
    big_mesh.topology.create_connectivity(cdim-1,cdim)
    bfacets_nodes_in_parent = indices_to_function(v_big,
                                                  subfacet_map[bfacets],
                                                  1,
                                                  name="bfacets_subproblem")
    partition = fem.Function(dg0_big,name="rank")
    partition.x.array[:] = rank
    writer = io.VTKFile(big_mesh.comm, f"big_mesh.pvd", "wb")
    writer.write_function([bfacets_nodes_in_parent,partition])
    writer.close()

def indices_to_function(space, indices, dim, name="f", remote=True, f = None):
    dofs = fem.locate_dofs_topological(space, dim, indices,remote)
    if f is None:
        f = fem.Function(space,name=name)
    else:
        f.x.array[:] = 0
    f.x.array[dofs] = 1
    return f

if __name__ ==" __main__":
    main()

```

1 processor (expected result):  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/9/9f9a2f102dd6360d4748219dd1ea36e0482e8a9b.png)

2 processors:  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/0/04cb5618c00d7a0a62c6b346b6341927f9bd503e.png)

I’m creating a sub-index map from the parent mesh and setting it as the facet index map for the submesh. I get the boundary facets from the submesh and I then try to get those boundary facets in the parent mesh using the facet correspondence provided by `create_subindex_map`.

I’m using the nightly docker image.

Any help is appreciated.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [July 1, 2024, 3:51pm UTC](https://fenicsproject.discourse.group/t/facet-map-from-submesh-to-parent-mesh/15118/2 "2024-07-01T15:51:53Z")

</div>

Is there a specific reason for wanting this map?  
If you want maps between parent and sub mesh indices, i would suggest having a look at:

> <https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py-L124>

as it shows how to map from parent to submesh entities. One could of course invert this map.

---

<div class="post-metadata">

**Author:** ![ordinary-slim](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ordinary-slim/32/4821_2.png) [@ordinary-slim](https://fenicsproject.discourse.group/u/ordinary-slim)\
**Post date:** [July 2, 2024, 10:25am UTC](https://fenicsproject.discourse.group/t/facet-map-from-submesh-to-parent-mesh/15118/3 "2024-07-02T10:25:43Z")

</div>

Thank you, your method did work. Here is the updated MVP:

```auto
from mpi4py import MPI
from dolfinx import mesh, fem, cpp, io
import numpy as np

comm = MPI.COMM_WORLD
rank = comm.Get_rank()

def build_subentity_to_parent_mapping(edim:int,pdomain:mesh.Mesh,cdomain:mesh.Mesh,subcell_map,subvertex_map):
    tdim = cdomain.topology.dim
    centity_map = cdomain.topology.index_map(edim)
    num_cents = centity_map.size_local + centity_map.num_ghosts
    subentity_map = np.full(num_cents,-1,dtype=np.int32)

    cdomain.topology.create_connectivity(edim,tdim)
    ccon_e2c = cdomain.topology.connectivity(edim,tdim)
    ccon_e2v = cdomain.topology.connectivity(edim,0)

    pdomain.topology.create_connectivity(tdim,edim)
    pcon_e2v = pdomain.topology.connectivity(edim,0)
    pcon_c2e = pdomain.topology.connectivity(tdim,edim)

    for ient in range(num_cents):
        entity_found = False
        cicells = ccon_e2c.links(ient)#child incident cells
        picells = subcell_map[cicells]#parent incident cells
        pinodes = subvertex_map[ccon_e2v.links(ient)]#parent incident nodes
        pinodes.sort()
        for picell in picells:
            if not(entity_found):
                for pient in pcon_c2e.links(picell):
                    pinodes2compare = pcon_e2v.links(pient).copy()
                    pinodes2compare.sort()
                    entity_found = (pinodes==pinodes2compare).all()
                    if entity_found:
                        subentity_map[ient] = pient
                        break
    return subentity_map

def main():
    points_side = 16
    big_mesh = mesh.create_unit_square(MPI.COMM_WORLD,
                                        points_side,
                                        points_side,
                                        )
    cdim = big_mesh.topology.dim
    v_big = fem.functionspace(big_mesh, ("CG", 1),)
    dg0_big = fem.functionspace(big_mesh, ("Discontinuous Lagrange", 0),)
    left_els = fem.locate_dofs_geometrical(dg0_big, lambda x : x[0] <= 0.5 )
    submesh, subcell_map, subvertex_map, _ = mesh.create_submesh(big_mesh,cdim,left_els)

    submesh.topology.create_entities(cdim-1)
    subfacet_map = build_subentity_to_parent_mapping(cdim-1,big_mesh,submesh,subcell_map,subvertex_map)

    bfacets = mesh.locate_entities_boundary(submesh,1,lambda x : np.ones(x.shape[1]))
    big_mesh.topology.create_connectivity(cdim-1,cdim)
    bfacets_nodes_in_parent = indices_to_function(v_big,
                                                  subfacet_map[bfacets],
                                                  1,
                                                  name="bfacets_subproblem")
    partition = fem.Function(dg0_big,name="rank")
    partition.x.array[:] = rank
    writer = io.VTKFile(big_mesh.comm, f"big_mesh.pvd", "wb")
    writer.write_function([bfacets_nodes_in_parent,partition])
    writer.close()

def indices_to_function(space, indices, dim, name="f", remote=True, f = None):
    dofs = fem.locate_dofs_topological(space, dim, indices,remote)
    if f is None:
        f = fem.Function(space,name=name)
    else:
        f.x.array[:] = 0
    f.x.array[dofs] = 1
    return f

if __name__ ==" __main__":
    main()

```