Coupling of two subdomains with new mixed-dimensional assembly

dolfinx
1

I am using the stable docker container for v0.9.0.:
DOLFINx version: 0.9.0 based on GIT commit: 4116cca8cf91f1a7e877d38039519349b3e08620

Hello,

I try to use the new mixed-dimensional assembly to couple two subdomains and run into a problem.
Maybe someone can clarify if I make a mistake or if the functionality I want is not (yet) implemented.

Consider a square domain cut in the middle (conforming with the mesh). I want to solve different equations on both subdomains and couple them in the interface weakly. Consider the following MWE:

def create_box(h):
    gmsh.initialize()
    mesh_comm = MPI.COMM_WORLD
    model_rank = 0
    if mesh_comm.rank == model_rank:
        gmsh.model.occ.addPoint(0, 0, 0, 1, 1)
        gmsh.model.occ.addPoint(0.5, 0, 0, 1, 2)
        gmsh.model.occ.addPoint(1, 0, 0, 1, 3)
        gmsh.model.occ.addPoint(1, 1, 0, 1, 4)
        gmsh.model.occ.addPoint(0.5, 1, 0, 1, 5)
        gmsh.model.occ.addPoint(0, 1, 0, 1, 6)
        gmsh.model.occ.addLine(1, 2, 1)
        gmsh.model.occ.addLine(2, 3, 2)
        gmsh.model.occ.addLine(3, 4, 3)
        gmsh.model.occ.addLine(4, 5, 4)
        gmsh.model.occ.addLine(5, 6, 5)
        gmsh.model.occ.addLine(6, 1, 6)
        gmsh.model.occ.addLine(2, 5, 7)
        gmsh.model.occ.addCurveLoop([1, 7, 5, 6], 1)
        gmsh.model.occ.addCurveLoop([2, 3, 4, 7], 2)
        gmsh.model.occ.addPlaneSurface([1], 1)
        gmsh.model.occ.addPlaneSurface([2], 2)
        gmsh.model.occ.synchronize()
        gmsh.model.addPhysicalGroup(dim=1, tags=[7], tag=1, name="interface")
        gmsh.model.addPhysicalGroup(dim=2, tags=[1], tag=1, name="left")
        gmsh.model.addPhysicalGroup(dim=2, tags=[2], tag=2, name="right")
        gmsh.model.mesh.setOrder(1)
        gmsh.option.setNumber("Mesh.MeshSizeMin", h)
        gmsh.option.setNumber("Mesh.MeshSizeMax", h)
        gmsh.model.mesh.generate(2)
    domain, cell_markers, facet_markers = io.gmshio.model_to_mesh(gmsh.model,
                                                               mesh_comm,
                                                               model_rank,
                                                               gdim=2)
    gmsh.finalize()
    return domain, cell_markers, facet_markers

This is the code I use to create my very simple domain. A square which is divided from middle top to bottom. This line is marked with facet_marker 1, the left subdomain with cell_marker 1 and the right subdomain with cell_marker 2.

Given the domain function, I run the following:

domain, cell_markers, facet_markers = create_box(0.5)
ddim = domain.geometry.dim
cell_name = domain.topology.cell_name()

P1 = basix.ufl.element("Lagrange", cell_name, 1)

# Define the subdomains
domain_1, entity_map_1, vertex_map_1, node_map_1 = mesh.create_submesh(domain, ddim, cell_markers.find(1))
domain_2, entity_map_2, vertex_map_2, node_map_2 = mesh.create_submesh(domain, ddim, cell_markers.find(2))

# Entity Maps
cell_imap_1 = domain.topology.index_map(2)
num_cells_1 = cell_imap_1.size_local + cell_imap_1.num_ghosts
parent_to_sub_1 = np.full(num_cells_1, -1)
parent_to_sub_1[entity_map_1] = np.arange(len(entity_map_1))

cell_imap_2 = domain.topology.index_map(2)
num_cells_2 = cell_imap_2.size_local + cell_imap_2.num_ghosts
parent_to_sub_2 = np.full(num_cells_2, -1)
parent_to_sub_2[entity_map_2] = np.arange(len(entity_map_2))

entity_maps = {domain_1:parent_to_sub_1, domain_2:parent_to_sub_2}

# Function Spaces
V = fem.functionspace(domain, P1)
V_1 = fem.functionspace(domain_1, P1)
V_2 = fem.functionspace(domain_2, P1)

dx = ufl.Measure("cell", domain, subdomain_data=cell_markers)
# Careful! Potentially inconsistent normal direction on interior facets!
# Just for illustration purposes
dS = ufl.Measure("interior_facet", domain,
                    subdomain_data=facet_markers)
n = ufl.FacetNormal(domain)

V_1 = fem.functionspace(domain_1, P1)
V_2 = fem.functionspace(domain_2, P1)
V = ufl.MixedFunctionSpace(V_1, V_2)

u1, u2 = ufl.TrialFunctions(V)
v1, v2 = ufl.TestFunctions(V)

# Bilinear form
a = u1 * v1 * dx(1) + u2 * v2 * dx(2)

# Comment this line for a working form
a += ufl.dot(ufl.grad(u1('+')), n('+')) * v2('-') * dS(1)

L = (fem.Constant(domain, PETSc.ScalarType(5)) * v1 * dx(1)
     + fem.Constant(domain, PETSc.ScalarType(0)) * v2 * dx(2))

a_blocked = ufl.extract_blocks(a)
L_blocked = ufl.extract_blocks(L)

a_blocked_compiled = fem.form(a_blocked, entity_maps=entity_maps)
L_blocked_compiled = fem.form(L_blocked, entity_maps=entity_maps)

A = fem.petsc.create_matrix_nest(a_blocked_compiled)

I am aware that my interior facet markers are of inconsistent normal direction, but that is not the point here. When I add an integral over the interface (dS(1)) that is coupling the two domains (test function from one domain, trial function from another), I get the PETSc error message
RuntimeError: Cannot insert rows that do not exist in the IndexMap. when creating the nested or blocked matrix.
When I comment this part, everything is fine.
(the same error occurs when I instead add a += u1 * v2 * dx to the bilinear form, not really physical, but omitting the interior facet integration)

This is closely related to this very nice tutorial, where the author also couples trial and test functions from different domains (in the nonlinear setting, so essentialy fem.Functions and ufl.TestFunctions), what seems to work just fine. The difference is that he couples an exterior subdomain (dim 1) to a domain (dim 2), and I need to couple a subdomain (dim 2) to another subdomain (dim 2) via the shared interface.

What am I overlooking? (My current implementation based on test and trial functions defined on the full domain with Dirichlet BCs to remove unwanted DOFs is working just fine.)

If I need to clarify something, please ask. Thank you very much!

2

I would consider

and the demo I made in

which couples two poisson problems discretised with different finite element spaces over a common boundary.

I have a similar implementation for RTC at

1 Like
3

Thank you for your suggestion. The problem is indeed solved in your examples. For future reference:

In the code, the entity_maps have to be adjusted like this (after creation):

domain.topology.create_connectivity(ddim-1, ddim)
facet_to_cell = domain.topology.connectivity(ddim-1, ddim)
for facet in facet_markers.find(1):
    cells = facet_to_cell.links(facet)
    map_1 = parent_to_sub_1[cells]
    map_2 = parent_to_sub_2[cells]
    parent_to_sub_1[cells] = max(map_1)
    parent_to_sub_2[cells] = max(map_2)

Otherwise, the interior facet integration with the global dS measure is not defined properly (because the global interior facets become exterior facets on the submesh). The modification of the entity_maps adds the corresponding cells from the other submesh to the entity_maps, such that integration is then well-defined.