Integrating over interface between two functionspaces

dolfinx
1

I want to create two functionspaces from one mesh, e.g. as in this image, with the interface \Gamma between the subdomains \Omega_A and \Omega_B at x=0.

Bildschirmfoto 2024-06-20 um 15.42.21
Bildschirmfoto 2024-06-20 um 15.42.21732×526 103 KB

I’ve used create_submesh to create two seperate meshes for each domain. Assembling vectors inside the domain works. Now I want to assemble Neumann BCs that are functions of the unknowns on both sides of the interface, e.g. f = 0.5(u_B - u_a). Think chemical reactions between species in \Omega_A and \Omega_B.

I created a Measure ds with the marker 1 at the interface, and tried to assemble \int_{\Gamma} v_a f \mathrm{d}s. This fails with the following error message:

IndexError                                Traceback (most recent call last)
File /home/ubuntu/interface_example.py:72
     68 boundary_term = v_a * f * ds(1)
     70 rhs = form(boundary_term)
---> 72 rhs_vector = assemble_vector(rhs)

File /usr/lib/python3.12/functools.py:909, in singledispatch.<locals>.wrapper(*args, **kw)
    905 if not args:
    906     raise TypeError(f'{funcname} requires at least '
    907                     '1 positional argument')
--> 909 return dispatch(args[0].__class__)(*args, **kw)

File /usr/local/dolfinx-real/lib/python3.12/dist-packages/dolfinx/fem/assemble.py:180, in _assemble_vector_form(L, constants, coeffs)
    178 b.array[:] = 0
    179 constants = constants or _pack_constants(L._cpp_object)
--> 180 coeffs = coeffs or _pack_coefficients(L._cpp_object)
    181 _assemble_vector_array(b.array, L, constants, coeffs)
    182 return b

IndexError: map::at

A M(n)WE of all I have described so far:

#!/usr/bin/env python3
from mpi4py import MPI
import numpy as np

import gmsh

from ufl import Measure, TestFunction
from dolfinx.io.gmshio import model_to_mesh
from dolfinx.mesh import create_submesh, locate_entities, meshtags
from dolfinx.fem import functionspace, Function, form, assemble_vector

def create_model():
    gmsh.initialize()
    gmsh.model.add("t2")

    lc = 0.25
    gmsh.model.occ.addPoint(1, 0, 0, lc, 1)
    gmsh.model.occ.addPoint(1, 1, 0, lc, 2)
    gmsh.model.occ.addPoint(1, 1, 1, lc, 3)
    gmsh.model.occ.addPoint(1, 0, 1, lc, 4)
    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, 1, 4)
    gmsh.model.occ.addCurveLoop([1, 2, 3, 4], 1)
    gmsh.model.occ.addPlaneSurface([1], 1)

    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(3, [1], name="left_volume")
    gmsh.model.addPhysicalGroup(3, [2], name="right_volume")

    gmsh.model.mesh.generate(3)

    return gmsh.model

mesh, cell_tags, facet_tags = model_to_mesh(create_model(), MPI.COMM_WORLD, rank=0)

left_indices = np.where(cell_tags.values==1)[0]
right_indices = np.where(cell_tags.values==2)[0]

submesh_left, cell_map_left, _, _ = create_submesh(mesh, 3, left_indices)
submesh_right, cell_map_right, _, _ = create_submesh(mesh, 3, right_indices)

V_A = functionspace(submesh_left, ("Lagrange", 1))
V_B = functionspace(submesh_right, ("Lagrange", 1))

u_a = Function(V_A)
u_b = Function(V_B)
v_a = TestFunction(V_A)
v_b = TestFunction(V_B)

u_a.vector.array.fill(1.0)
u_b.vector.array.fill(2.0)

interface_facets = locate_entities(submesh_left, 2, lambda x: np.isclose(x[0], 0.0))
marker = np.full_like(interface_facets, 1)
facet_tags = meshtags(submesh_left, 2, interface_facets, marker)
ds = Measure("ds", domain=submesh_left, subdomain_data=facet_tags)

f = 0.5*(u_b - u_a)
boundary_term = v_a * f * ds(1)

rhs = form(boundary_term)

rhs_vector = assemble_vector(rhs)

This leads me to the following questions:

  1. Is this in general the correct approach?
  2. I found that I can supply form() with entity_maps, which I think are the maps this index error occurs, but I’m not sure what to put in there. The submeshes don’t have any cells in common; would supplying facet and point maps suffice? How should they be constructed?
  3. Can I put this in more general terms, so that I don’t have to do the interface tags & measure twice? For a form v_b * f * ds(1) I would presumeably need a different ds, right?

I’m running MacOS 14 on ARM, and tried with both Dolfinx 0.8 from conda, and nightly using docker.

2

This exact thing has been covered in: Poisson submesh DG coupling · GitHub
where I do this for a 2D problem (with coupling terms, nitsche/DG style) across the interface.

One needs to use dS integrals (a bit clunky I would admit) to get the mappings correctly.
Another place to look is:
mixed_domain_demos/poisson_domain_decomp.py at main · jpdean/mixed_domain_demos · GitHub by @jpdean

3

I will have a look, thank you.

4

Thanks again for the helpful examples. I have a couple of follow up questions regarding your problem.py example you linked to in the gist:

  1. In lines 238 & 240, where you are creating the “inverse” entity maps, shouldn’t they be of size num_cells_local, rather than num_facets_local?
  2. I don’t quite understand the restrictions (lines 304ff). I’m assuming ("+") means the side of the facet the normal vector is pointing.
    2a) These are continuous spaces, so what difference does it make?
    2b) In any case, I would have assumed it is relative to each space, so why not “outside” ("+")for both? Or relative to the interface orientation, in which case how do I find out in which direction the interface normal is pointing?
5

Yes, that is true. It does work atm as there are more facets than cells on the process.
I’ll change that.

This is to work around the fact that one has to use the dS measure to define the mix terms.
As each submesh only has a common interface, no shared cells, we need to pass information about cells in each submesh to the integration kernel.
This is supported through the dS measure.
As we would like to have terms from both sides at the same time, we choose one to be +, the other minus. This makes sure that we get all contributions along the common boundary. It doesn’t matter which side of the interface gets + or -.

6

Thank you, that was very helpful.