Facet map from submesh to parent mesh

dolfinx
1

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:

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

2 processors:
image

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.

2

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:

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

1 Like
3

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

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()