DBCs on submesh yield to occasional segfault

General
1

Hello,

Consider the following MWE that is producing a segfault in latest dolfinx nightly (pulled 17 Jun) occasionally when trying to set DBCs on a submesh (though not on every call, hence the locate_dofs_topological is called in a loop of 100 times):

#!/usr/bin/env python3
from mpi4py import MPI
import numpy as np
from dolfinx import io, fem, mesh
import ufl

comm = MPI.COMM_WORLD

dbcs_on_submeshes = True # True, False

mshfile_d='artseg-fsi-quad_domain.xdmf'
mshfile_b='artseg-fsi-quad_boundary.xdmf'

with io.XDMFFile(comm, mshfile_d, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    msh = infile.read_mesh(name="Grid")
    mt_d = infile.read_meshtags(msh, name="Grid")

msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
with io.XDMFFile(comm, mshfile_b, 'r', encoding=io.XDMFFile.Encoding.ASCII) as infile:
    mt_b1 = infile.read_meshtags(msh, name="Grid")

# cf. https://fenicsproject.discourse.group/t/transfer-meshtags-to-submesh-in-dolfinx/8952/6
def meshtags_parent_to_child(mshtags, childmsh, childmsh_emap, parentmsh):
    dim = parentmsh.topology.dim-1
    d_map = parentmsh.topology.index_map(dim)
    all_ent = d_map.size_local + d_map.num_ghosts
    # create array with zeros for all entities that are not marked
    all_values = np.zeros(all_ent, dtype=np.int32)
    all_values[mshtags.indices] = mshtags.values
    c_to_e = parentmsh.topology.connectivity(parentmsh.topology.dim, dim)
    childmsh.topology.create_entities(dim)
    subf_map = childmsh.topology.index_map(dim)
    childmsh.topology.create_connectivity(parentmsh.topology.dim, dim)
    c_to_e_sub = childmsh.topology.connectivity(parentmsh.topology.dim, dim)
    num_sub_ent = subf_map.size_local + subf_map.size_global
    sub_values = np.empty(num_sub_ent, dtype=np.int32)
    for i, entity in enumerate(childmsh_emap):
        parent_ent = c_to_e.links(entity)
        child_ent = c_to_e_sub.links(i)
        for child, parent in zip(child_ent, parent_ent):
            sub_values[child] = all_values[parent]
    return mesh.meshtags(childmsh, childmsh.topology.dim-1, np.arange(num_sub_ent, dtype=np.int32), sub_values)

for i in range(100):

    msh_emap_solid = mesh.create_submesh(msh, msh.topology.dim, mt_d.indices[mt_d.values == 1])[0:2]
    msh_emap_fluid = mesh.create_submesh(msh, msh.topology.dim, mt_d.indices[mt_d.values == 2])[0:2]

    mt_b1_fluid = meshtags_parent_to_child(mt_b1, msh_emap_fluid[0], msh_emap_fluid[1], msh)
    mt_b1_solid = meshtags_parent_to_child(mt_b1, msh_emap_solid[0], msh_emap_solid[1], msh)

    # fluid and solid function spaces
    P_v = ufl.VectorElement("CG", msh_emap_fluid[0].ufl_cell(), 2)
    V_v = fem.FunctionSpace(msh_emap_fluid[0], P_v)
    P_u = ufl.VectorElement("CG", msh_emap_solid[0].ufl_cell(), 2)
    V_u = fem.FunctionSpace(msh_emap_solid[0], P_u)

    if dbcs_on_submeshes:
        # fluid DBCs
        func_f = fem.Function(V_v) # zero function
        dbcs_fluid = []
        dbcs_fluid.append( fem.dirichletbc(func_f, fem.locate_dofs_topological(V_v, msh_emap_fluid[0].topology.dim-1, mt_b1_fluid.indices[mt_b1_fluid.values == 4])) )

        dofs_x_f = fem.locate_dofs_topological(V_v.sub(0), msh_emap_fluid[0].topology.dim-1, mt_b1_fluid.indices[mt_b1_fluid.values == 5])
        dbcs_fluid.append( fem.dirichletbc(func_f.sub(0), dofs_x_f) )

        # solid DBCs
        func_s = fem.Function(V_u) # zero function
        dbcs_solid = []
        dbcs_solid.append( fem.dirichletbc(func_s, fem.locate_dofs_topological(V_u, msh_emap_solid[0].topology.dim-1, mt_b1_solid.indices[mt_b1_solid.values == 4])) )

        dofs_x_s = fem.locate_dofs_topological(V_u.sub(0), msh_emap_solid[0].topology.dim-1, mt_b1_solid.indices[mt_b1_solid.values == 5])
        dbcs_solid.append( fem.dirichletbc(func_s.sub(0), dofs_x_s) )

    # some function spaces on full mesh
    P_ = ufl.VectorElement("CG", msh.ufl_cell(), 2)
    V_ = fem.FunctionSpace(msh, P_v)
    func_ = fem.Function(V_) # zero function
    
    dbcs_ = []
    dbcs_.append( fem.dirichletbc(func_, fem.locate_dofs_topological(V_, msh.topology.dim-1, mt_b1.indices[mt_b1.values == 5])) )
    dofs_x_ = fem.locate_dofs_topological(V_.sub(0), msh.topology.dim-1, mt_b1.indices[mt_b1.values == 5])
    dbcs_.append( fem.dirichletbc(func_.sub(0), dofs_x_) )


## check that mesh tags are actually written correctly
#with io.XDMFFile(comm, "sub_tag_solid.xdmf", "w") as xdmf:
    #xdmf.write_mesh(msh_emap_solid[0])
    #msh_emap_solid[0].topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
    #xdmf.write_meshtags(mt_b1_solid, msh_emap_solid[0].geometry)

#with io.XDMFFile(comm, "sub_tag_fluid.xdmf", "w") as xdmf:
    #xdmf.write_mesh(msh_emap_fluid[0])
    #msh_emap_fluid[0].topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
    #xdmf.write_meshtags(mt_b1_fluid, msh_emap_fluid[0].geometry)

Mesh files (hexa 27) can be retrieved from
http://www.deepambit.com/stuff/artseg-fsi-quad_domain.xdmf
http://www.deepambit.com/stuff/artseg-fsi-quad_boundary.xdmf

Error says
python3: /src/dolfinx/cpp/dolfinx/graph/AdjacencyList.h:102: int dolfinx::graph::AdjacencyList<T>::num_links(int) const [with T = int]: Assertion (node + 1) < (int)_offsets.size()’ failed.`

Setting dbcs_on_submeshes = False, hence bypassing that step and setting the DBCs on the main mesh works. Uncomment the last lines to verify that the meshtags are correctly transformed.

Am I doing something wrong here?

Thanks for any help.

Best,
Marc

2

I made a few changes to your code (simplifying the loop a bit):

# cf. https://fenicsproject.discourse.group/t/transfer-meshtags-to-submesh-in-dolfinx/8952/6
def meshtags_parent_to_child(mshtags, childmsh, childmsh_emap, parentmsh):
    dim = parentmsh.topology.dim-1
    d_map = parentmsh.topology.index_map(dim)
    all_ent = d_map.size_local + d_map.num_ghosts
    # create array with zeros for all entities that are not marked
    all_values = np.zeros(all_ent, dtype=np.int32)
    all_values[mshtags.indices] = mshtags.values
    c_to_e = parentmsh.topology.connectivity(parentmsh.topology.dim, dim)
    childmsh.topology.create_entities(dim)
    subf_map = childmsh.topology.index_map(dim)
    childmsh.topology.create_connectivity(parentmsh.topology.dim, dim)
    c_to_e_sub = childmsh.topology.connectivity(parentmsh.topology.dim, dim)
    num_sub_ent = subf_map.size_local + subf_map.size_global
    sub_values = np.zeros(num_sub_ent, dtype=np.int32)
    for i, entity in enumerate(childmsh_emap):
        parent_ent = c_to_e.links(entity)
        child_ent = c_to_e_sub.links(i)
        for child, parent in zip(child_ent, parent_ent):
            sub_values[child] = all_values[parent]
    print(num_sub_ent, sub_values)
    return mesh.meshtags(childmsh, childmsh.topology.dim-1, np.arange(num_sub_ent, dtype=np.int32), sub_values)

for i in range(100):
    msh_emap_fluid = mesh.create_submesh(msh, msh.topology.dim, mt_d.indices[mt_d.values == 2])[0:2]
    mt_b1_fluid = meshtags_parent_to_child(mt_b1, msh_emap_fluid[0], msh_emap_fluid[1], msh)

    # fluid and solid function spaces
    P_v = ufl.VectorElement("CG", msh_emap_fluid[0].ufl_cell(), 2)
    V_v = fem.FunctionSpace(msh_emap_fluid[0], P_v)

    if dbcs_on_submeshes:
        # fluid DBCs
        func_f = fem.Function(V_v) # zero function
        dbcs_fluid = []
        dbcs_fluid.append( fem.dirichletbc(func_f, fem.locate_dofs_topological(V_v, msh_emap_fluid[0].topology.dim-1, mt_b1_fluid.indices[mt_b1_fluid.values == 4])) )
        V0_v, _ = V_v.sub(0).collapse()
        func_f_x = fem.Function(V0_v)
        dofs_x_f = fem.locate_dofs_topological((V_v.sub(0), V0_v), msh_emap_fluid[0].topology.dim-1, mt_b1_fluid.indices[mt_b1_fluid.values == 5])
        dbcs_fluid.append( fem.dirichletbc(func_f_x, dofs_x_f, V_v.sub(0)) )

Most noteably, I changed

to
sub_values = np.zeros(num_sub_ent, dtype=np.int32)
and

to


        V0_v, _ = V_v.sub(0).collapse()
        func_f_x = fem.Function(V0_v)
        print( mt_b1_fluid.indices[mt_b1_fluid.values == 5])
        dofs_x_f = fem.locate_dofs_topological((V_v.sub(0), V0_v), msh_emap_fluid[0].topology.dim-1, mt_b1_fluid.indices[mt_b1_fluid.values == 5])
        dbcs_fluid.append( fem.dirichletbc(func_f_x, dofs_x_f, V_v.sub(0)) )
3

Also, with similar changes to solid, the code runs reliably:


        # solid DBCs
        func_s = fem.Function(V_u) # zero function
        dbcs_solid = []
        dbcs_solid.append( fem.dirichletbc(func_s, fem.locate_dofs_topological(V_u, msh_emap_solid[0].topology.dim-1, mt_b1_solid.indices[mt_b1_solid.values == 4])) )
        V0_u, _ = V_u.sub(0).collapse()
        func_u_x = fem.Function(V0_u)
        dofs_x_s = fem.locate_dofs_topological((V_u.sub(0), V0_u), msh_emap_solid[0].topology.dim-1, mt_b1_solid.indices[mt_b1_solid.values == 5])
        dbcs_solid.append( fem.dirichletbc(func_u_x, dofs_x_s, V_u.sub(0)) )
4

Awesome, thanks, that worked!