# Integrating over interface between two functionspaces **URL:** \ **Category:** dolfinx\ **Created:** [June 20, 2024, 2:32pm UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037 "2024-06-20T14:32:00Z")\ **Posts on this page:** 12\ **Page:** 1
**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\ **Post date:** [June 20, 2024, 2:32pm UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/1 "2024-06-20T14:32:00Z")
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](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6a7768609b3a6d3e0c16eb61b170ad24df739ab4.png) 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: ```auto 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..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: ```python #!/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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [June 20, 2024, 3:07pm UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/2 "2024-06-20T15:07:49Z")
This exact thing has been covered in: [Poisson submesh DG coupling · GitHub](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py) 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](https://github.com/jpdean/mixed_domain_demos/blob/main/poisson_domain_decomp.py) by @jpdean ---
**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\ **Post date:** [June 20, 2024, 3:33pm UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/3 "2024-06-20T15:33:00Z")
I will have a look, thank you. ---
**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\ **Post date:** [June 24, 2024, 9:44am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/4 "2024-06-24T09:44:27Z")
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? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [June 24, 2024, 6:14pm UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/5 "2024-06-24T18:14:04Z")
> [@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 -. ---
**Author:** ![Psirus](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/psirus/32/624_2.png) [@Psirus](https://fenicsproject.discourse.group/u/Psirus)\ **Post date:** [June 25, 2024, 7:07am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/6 "2024-06-25T07:07:55Z")
Thank you, that was very helpful. ---
**Author:** ![Shi](https://avatars.discourse-cdn.com/v4/letter/s/f14d63/32.png) [@Shi](https://fenicsproject.discourse.group/u/Shi)\ **Post date:** [December 25, 2024, 7:16am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/7 "2024-12-25T07:16:56Z")
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? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [December 25, 2024, 7:31am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/8 "2024-12-25T07:31:43Z")
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](https://github.com/scientificcomputing/scifem/blob/main/src/scifem/solvers.py#L248) and usage in [scifem/tests/test\_blocked\_newton\_solver.py at main · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/blob/main/tests/test_blocked_newton_solver.py#L70) ---
**Author:** ![Shi](https://avatars.discourse-cdn.com/v4/letter/s/f14d63/32.png) [@Shi](https://fenicsproject.discourse.group/u/Shi)\ **Post date:** [December 25, 2024, 8:46am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/9 "2024-12-25T08:46:18Z")
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](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py) or the BlockedNewtonSolver in scifem that you provided [scifem/src/scifem/solvers.py at main · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/blob/main/src/scifem/solvers.py#L248). 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. ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [December 25, 2024, 8:59am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/10 "2024-12-25T08:59:39Z")
> [@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](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py) or the BlockedNewtonSolver in scifem that you provided [scifem/src/scifem/solvers.py at main · scientificcomputing/scifem · GitHub](https://github.com/scientificcomputing/scifem/blob/main/src/scifem/solvers.py#L248). 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. ---
**Author:** ![Shi](https://avatars.discourse-cdn.com/v4/letter/s/f14d63/32.png) [@Shi](https://fenicsproject.discourse.group/u/Shi)\ **Post date:** [December 25, 2024, 9:12am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/11 "2024-12-25T09:12:22Z")
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? ---
**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\ **Post date:** [December 25, 2024, 9:29am UTC](https://fenicsproject.discourse.group/t/integrating-over-interface-between-two-functionspaces/15037/12 "2024-12-25T09:29:22Z")
> [@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).