Integrating over interface between two functionspaces

[Psirus]

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.

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.

[dokken]

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

[Psirus]

I will have a look, thank you.

[Psirus]

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?

[dokken]

@Psirus

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?

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

@Psirus

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?

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 -.

[Psirus]

Thank you, that was very helpful.

[Shi]

hello ,After reading the code you wrote, I have a question. Is it necessary to rewrite the NewtonSolver to meet the requirement of solving on multiple submesh just like what do in that code? I’m not quite clear about this issue. Could you please help me solve it?

[dokken]

Yes,
There are two implementations of such a solver in scifem. I would recommend the blockednewtonsolver: scifem/src/scifem/solvers.py at main · scientificcomputing/scifem · GitHub
and usage in
scifem/tests/test_blocked_newton_solver.py at main · scientificcomputing/scifem · GitHub

[Shi]

Thank you for your reply. Currently, I have several questions regarding this discontinuous problem.

1. I’m using version 0.9 of Dolfinx. I’m wondering whether it would be better to use the solver you wrote in Poisson submesh DG coupling · GitHub or the BlockedNewtonSolver in scifem that you provided scifem/src/scifem/solvers.py at main · scientificcomputing/scifem · GitHub. Please forgive me as I’m a novice. I don’t quite understand why the traditional Newton solver is not very suitable for this kind of problem.
2. For the Dirichlet boundary conditions in problems with discontinuities at interfaces, should I apply them using the function space in the submesh or the function space in the parent mesh?

Thank you for your answers.

[dokken]

@Shi

I’m using version 0.9 of Dolfinx. I’m wondering whether it would be better to use the solver you wrote in Poisson submesh DG coupling · GitHub or the BlockedNewtonSolver in scifem that you provided scifem/src/scifem/solvers.py at main · scientificcomputing/scifem · GitHub. Please forgive me as I’m a novice. I don’t quite understand why the traditional Newton solver is not very suitable for this kind of problem.

The BlockedNewtonSolver is recommended.
It builds on top of the dolfinx Newton solver, extending the assembly routines to use assemble_matrix_block rather than assemble_matrix.

@Shi

For the Dirichlet boundary conditions in problems with discontinuities at interfaces, should I apply them using the function space in the submesh or the function space in the parent mesh?

As your unknown (trial function) lives on the submesh, the dirichletbc has to be associated with the function space on that mesh.

[Shi]

Thank you for your answer. In the mixed Poisson equation you dealt with before, when defining the boundary conditions, you first transformed the mesh tags from the parent mesh to the submesh, and then found the entities of the boundary marked by the mesh tags in the submesh. Is this the correct understanding?

Since my test and trial functions are all on the submesh, why do the integral measures ds and dx need to be defined on the parent mesh?

[dokken]

@Shi

Thank you for your answer. In the mixed Poisson equation you dealt with before, when defining the boundary conditions, you first transformed the mesh tags from the parent mesh to the submesh, and then found the entities of the boundary marked by the mesh tags in the submesh. Is this the correct understanding?

Yes.

@Shi

Since my test and trial functions are all on the submesh, why do the integral measures ds and dx need to be defined on the parent mesh?

That depends on the variational form.

If you need to work with two different submeshes omega0, omega1 sharing an interface, you need the integration measures to be specified on the parent mesh.

However, if you only work on the submesh, with dx or ds integrals, you can use measures on the submesh, as long as you provide the correct entity map, mapping any cell of the submesh to the parent mesh (if you have coefficients or constants defined on the parent mesh in your form).