# Discontinuity at interface using mixed domains

**URL:** <https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040>\
**Category:** variational formulation\
**Tags:** dolfinx\
**Created:** [June 20, 2024, 4:56pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040 "2024-06-20T16:56:25Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [June 20, 2024, 4:56pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/1 "2024-06-20T16:56:25Z")

</div>

Hi all,

Exactly like in [this topic](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-at-interface-using-dg-method/11780/18) I am trying to impose a discontinuity at an interface between two subdomains. This time having two subdomains with only CG elements and coupling them with a `dS` term.

@dokken published [this gist](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py) that I am trying to adapt in order to have:

\frac{u(-)}{K(-)} = \frac{u(+)}{K(+)}

\nabla u(-) - \nabla u(+) = 0

However I am struggling to adapt the interface terms to achieve the desired result.  
This is the DG formulation we had earlier:

```python
# Interface 
F += - dot(avg(grad(v)), n('-'))*(u('-')*(K1/K2-1))*dS(2)
F += alpha/avg(h)*dot(jump(v,n),n('-'))*(u('-')*(K1/K2-1))*dS(2)

# symmetry
F += - dot(avg(grad(v)), jump(u, n))*dS(2)
# coercivity
F += + alpha/avg(h)*dot(jump(v, n), jump(u, n))*dS(2)

```

---

<div class="post-metadata">

**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, 7:27pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/2 "2024-06-20T19:27:39Z")

</div>

I think this is the variational form that you would like:

```python
# SPDX-License-Identifier: MIT

from mpi4py import MPI
import dolfinx
import dolfinx.fem.petsc
import ufl
import numpy as np
from petsc4py import PETSc

class NewtonSolver:
    max_iterations: int
    bcs: list[dolfinx.fem.DirichletBC]
    A: PETSc.Mat
    b: PETSc.Vec
    J: dolfinx.fem.Form
    b: dolfinx.fem.Form
    dx: PETSc.Vec

    def __init__ (
        self,
        F: list[dolfinx.fem.form],
        J: list[list[dolfinx.fem.form]],
        w: list[dolfinx.fem.Function],
        bcs: list[dolfinx.fem.DirichletBC] | None = None,
        max_iterations: int = 5,
        petsc_options: dict[str, str | float | int | None] = None,
        problem_prefix="newton",
    ):
        self.max_iterations = max_iterations
        self.bcs = [] if bcs is None else bcs
        self.b = dolfinx.fem.petsc.create_vector_block(F)
        self.F = F
        self.J = J
        self.A = dolfinx.fem.petsc.create_matrix_block(J)
        self.dx = self.A.createVecLeft()
        self.w = w
        self.x = dolfinx.fem.petsc.create_vector_block(F)

        # Set PETSc options
        opts = PETSc.Options()
        if petsc_options is not None:
            for k, v in petsc_options.items():
                opts[k] = v

        # Define KSP solver
        self._solver = PETSc.KSP().create(self.b.getComm().tompi4py())
        self._solver.setOperators(self.A)
        self._solver.setFromOptions()

        # Set matrix and vector PETSc options
        self.A.setFromOptions()
        self.b.setFromOptions()

    def solve(self, tol=1e-6, beta=1.0):
        i = 0

        while i < self.max_iterations:
            dolfinx.cpp.la.petsc.scatter_local_vectors(
                self.x,
                [si.x.petsc_vec.array_r for si in self.w],
                [
                    (
                        si.function_space.dofmap.index_map,
                        si.function_space.dofmap.index_map_bs,
                    )
                    for si in self.w
                ],
            )
            self.x.ghostUpdate(
                addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD
            )

            # Assemble F(u_{i-1}) - J(u_D - u_{i-1}) and set du|_bc= u_D - u_{i-1}
            with self.b.localForm() as b_local:
                b_local.set(0.0)
            dolfinx.fem.petsc.assemble_vector_block(
                self.b, self.F, self.J, bcs=self.bcs, x0=self.x, scale=-1.0
            )
            self.b.ghostUpdate(
                PETSc.InsertMode.INSERT_VALUES, PETSc.ScatterMode.FORWARD
            )

            # Assemble Jacobian
            self.A.zeroEntries()
            dolfinx.fem.petsc.assemble_matrix_block(self.A, self.J, bcs=self.bcs)
            self.A.assemble()

            self._solver.solve(self.b, self.dx)
            # self._solver.view()
            assert (
                self._solver.getConvergedReason() > 0
            ), "Linear solver did not converge"
            offset_start = 0
            for s in self.w:
                num_sub_dofs = (
                    s.function_space.dofmap.index_map.size_local
                    * s.function_space.dofmap.index_map_bs
                )
                s.x.petsc_vec.array_w[:num_sub_dofs] -= (
                    beta * self.dx.array_r[offset_start : offset_start + num_sub_dofs]
                )
                s.x.petsc_vec.ghostUpdate(
                    addv=PETSc.InsertMode.INSERT, mode=PETSc.ScatterMode.FORWARD
                )
                offset_start += num_sub_dofs
            # Compute norm of update

            correction_norm = self.dx.norm(0)
            print(f"Iteration {i}: Correction norm {correction_norm}")
            if correction_norm < tol:
                break
            i += 1

    def __del__ (self):
        self.A.destroy()
        self.b.destroy()
        self.dx.destroy()
        self._solver.destroy()
        self.x.destroy()

def transfer_meshtags_to_submesh(
    mesh, entity_tag, submesh, sub_vertex_to_parent, sub_cell_to_parent
):
    """
    Transfer a meshtag from a parent mesh to a sub-mesh.
    """

    tdim = mesh.topology.dim
    cell_imap = mesh.topology.index_map(tdim)
    num_cells = cell_imap.size_local + cell_imap.num_ghosts
    mesh_to_submesh = np.full(num_cells, -1)
    mesh_to_submesh[sub_cell_to_parent] = np.arange(
        len(sub_cell_to_parent), dtype=np.int32
    )
    sub_vertex_to_parent = np.asarray(sub_vertex_to_parent)

    submesh.topology.create_connectivity(entity_tag.dim, 0)

    num_child_entities = (
        submesh.topology.index_map(entity_tag.dim).size_local
        + submesh.topology.index_map(entity_tag.dim).num_ghosts
    )
    submesh.topology.create_connectivity(submesh.topology.dim, entity_tag.dim)

    c_c_to_e = submesh.topology.connectivity(submesh.topology.dim, entity_tag.dim)
    c_e_to_v = submesh.topology.connectivity(entity_tag.dim, 0)

    child_markers = np.full(num_child_entities, 0, dtype=np.int32)

    mesh.topology.create_connectivity(entity_tag.dim, 0)
    mesh.topology.create_connectivity(entity_tag.dim, mesh.topology.dim)
    p_f_to_v = mesh.topology.connectivity(entity_tag.dim, 0)
    p_f_to_c = mesh.topology.connectivity(entity_tag.dim, mesh.topology.dim)
    sub_to_parent_entity_map = np.full(num_child_entities, -1, dtype=np.int32)
    for facet, value in zip(entity_tag.indices, entity_tag.values):
        facet_found = False
        for cell in p_f_to_c.links(facet):
            if facet_found:
                break
            if (child_cell := mesh_to_submesh[cell]) != -1:
                for child_facet in c_c_to_e.links(child_cell):
                    child_vertices = c_e_to_v.links(child_facet)
                    child_vertices_as_parent = sub_vertex_to_parent[child_vertices]
                    is_facet = np.isin(
                        child_vertices_as_parent, p_f_to_v.links(facet)
                    ).all()
                    if is_facet:
                        child_markers[child_facet] = value
                        facet_found = True
                        sub_to_parent_entity_map[child_facet] = facet
    tags = dolfinx.mesh.meshtags(
        submesh,
        entity_tag.dim,
        np.arange(num_child_entities, dtype=np.int32),
        child_markers,
    )
    tags.name = entity_tag.name
    return tags, sub_to_parent_entity_map

def bottom_boundary(x):
    return np.isclose(x[1], 0.0)

def top_boundary(x):
    return np.isclose(x[1], 1.0)

def half(x):
    return x[1] <= 0.5 + 1e-14

mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 10, 10, dolfinx.mesh.CellType.triangle)

# Split domain in half and set an interface tag of 5
gdim = mesh.geometry.dim
tdim = mesh.topology.dim
fdim = tdim - 1
top_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, top_boundary)
bottom_facets = dolfinx.mesh.locate_entities_boundary(mesh, fdim, bottom_boundary)
num_facets_local = (
    mesh.topology.index_map(fdim).size_local + mesh.topology.index_map(fdim).num_ghosts
)
facets = np.arange(num_facets_local, dtype=np.int32)
values = np.full_like(facets, 0, dtype=np.int32)
values[top_facets] = 1
values[bottom_facets] = 2

bottom_cells = dolfinx.mesh.locate_entities(mesh, tdim, half)
num_cells_local = (
    mesh.topology.index_map(tdim).size_local + mesh.topology.index_map(tdim).num_ghosts
)
cells = np.full(num_cells_local, 4, dtype=np.int32)
cells[bottom_cells] = 3
ct = dolfinx.mesh.meshtags(
    mesh, tdim, np.arange(num_cells_local, dtype=np.int32), cells
)
all_b_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(3), tdim, fdim
)
all_t_facets = dolfinx.mesh.compute_incident_entities(
    mesh.topology, ct.find(4), tdim, fdim
)
interface = np.intersect1d(all_b_facets, all_t_facets)
values[interface] = 5

mt = dolfinx.mesh.meshtags(mesh, mesh.topology.dim - 1, facets, values)

submesh_b, submesh_b_to_mesh, b_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(3)
)[0:3]
submesh_t, submesh_t_to_mesh, t_v_map = dolfinx.mesh.create_submesh(
    mesh, tdim, ct.find(4)
)[0:3]
parent_to_sub_b = np.full(num_facets_local, -1, dtype=np.int32)
parent_to_sub_b[submesh_b_to_mesh] = np.arange(len(submesh_b_to_mesh), dtype=np.int32)
parent_to_sub_t = np.full(num_facets_local, -1, dtype=np.int32)
parent_to_sub_t[submesh_t_to_mesh] = np.arange(len(submesh_t_to_mesh), dtype=np.int32)    

# We need to modify the cell maps, as for `dS` integrals of interfaces between submeshes, there is no entity to map to.
# We use the entity on the same side to fix this (as all restrictions are one-sided)

# Transfer meshtags to submesh
ft_b, b_facet_to_parent = transfer_meshtags_to_submesh(
    mesh, mt, submesh_b, b_v_map, submesh_b_to_mesh
)
ft_t, t_facet_to_parent = transfer_meshtags_to_submesh(
    mesh, mt, submesh_t, t_v_map, submesh_t_to_mesh
)

t_parent_to_facet = np.full(num_facets_local, -1)
t_parent_to_facet[t_facet_to_parent] = np.arange(len(t_facet_to_parent), dtype=np.int32)

# Hack, as we use one-sided restrictions, pad dS integral with the same entity from the same cell on both sides
mesh.topology.create_connectivity(fdim, tdim)
f_to_c = mesh.topology.connectivity(fdim, tdim)
for facet in mt.find(5):
    cells = f_to_c.links(facet)
    assert len(cells) == 2
    b_map = parent_to_sub_b[cells]
    t_map = parent_to_sub_t[cells]
    parent_to_sub_b[cells] = max(b_map)
    parent_to_sub_t[cells] = max(t_map)

entity_maps = {submesh_b: parent_to_sub_b, submesh_t: parent_to_sub_t}
#entity_maps = {submesh_b._cpp_object: parent_to_sub_b, submesh_t._cpp_object: parent_to_sub_t}

def define_interior_eq(mesh,degree, submesh, submesh_to_mesh, value):
    # Compute map from parent entity to submesh cell
    codim = mesh.topology.dim - submesh.topology.dim
    ptdim = mesh.topology.dim - codim
    num_entities = (
        mesh.topology.index_map(ptdim).size_local
        + mesh.topology.index_map(ptdim).num_ghosts
    )
    mesh_to_submesh = np.full(num_entities, -1)
    mesh_to_submesh[submesh_to_mesh] = np.arange(len(submesh_to_mesh), dtype=np.int32)

    V = dolfinx.fem.functionspace(submesh, ("Lagrange", degree))
    u = dolfinx.fem.Function(V)
    v = ufl.TestFunction(V)
    ct_r = dolfinx.mesh.meshtags(mesh, mesh.topology.dim, submesh_to_mesh, np.full_like(submesh_to_mesh, 1, dtype=np.int32))
    val = dolfinx.fem.Constant(submesh, value)
    dx_r = ufl.Measure("dx", domain=mesh, subdomain_data=ct_r, subdomain_id=1)
    F = ufl.inner(ufl.grad(u), ufl.grad(v)) * dx_r - val * v * dx_r
    return u, F, mesh_to_submesh

u_0, F_00, m_to_b = define_interior_eq(mesh, 2, submesh_b, submesh_b_to_mesh, 0.0)
u_1, F_11, m_to_t = define_interior_eq(mesh, 1, submesh_t, submesh_t_to_mesh, 0.0)
u_0.name = "u_b"
u_1.name = "u_t"

# Add coupling term to the interface
# Get interface markers on submesh b
dInterface = ufl.Measure("dS", domain=mesh, subdomain_data=mt, subdomain_id=5)
b_res = "+"
t_res = "-"

v_b = ufl.TestFunction(u_0.function_space)(b_res)
v_t = ufl.TestFunction(u_1.function_space)(t_res)
u_b = u_0(b_res)
u_t = u_1(t_res)

def mixed_term(u, v, n):
    return ufl.dot(ufl.grad(u), n) * v

W_0 = dolfinx.fem.functionspace(submesh_b, ("DG", 0))
K_0 = dolfinx.fem.Function(W_0)
K_0.x.array[:] = 1
W_1 = dolfinx.fem.functionspace(submesh_t, ("DG", 0))
K_1 = dolfinx.fem.Function(W_1)
K_1.x.array[:] = 1.5

n = ufl.FacetNormal(mesh)
n_b = n(b_res)
n_t = n(t_res)
K_b = K_0(b_res)
K_t = K_1(t_res)
cr = ufl.Circumradius(mesh)
h_b = 2 * cr(b_res)
h_t = 2 * cr(t_res)
gamma = 50.0

F_0 = (
    -0.5 * mixed_term((u_b + u_t), v_b, n_b) * dInterface
    - 0.5 * mixed_term(v_b, (u_b/K_b - u_t/K_t), n_b) * dInterface
)

F_1 = (
    +0.5 * mixed_term((u_b + u_t), v_t, n_b) * dInterface
    - 0.5 * mixed_term(v_t, (u_b/K_b - u_t/K_t), n_b) * dInterface
)
F_0 += 2 * gamma / (h_b + h_t) * (u_b/K_b - u_t/K_t) * v_b * dInterface
F_1 += -2 * gamma / (h_b + h_t) * (u_b/K_b - u_t/K_t) * v_t * dInterface

F_0 += F_00
F_1 += F_11

jac00 = ufl.derivative(F_0, u_0)

jac01 = ufl.derivative(F_0, u_1)

jac10 = ufl.derivative(F_1, u_0)
jac11 = ufl.derivative(F_1, u_1)
J00 = dolfinx.fem.form(jac00, entity_maps=entity_maps)

J01 = dolfinx.fem.form(jac01, entity_maps=entity_maps)
J10 = dolfinx.fem.form(jac10, entity_maps=entity_maps)
J11 = dolfinx.fem.form(jac11, entity_maps=entity_maps)
J = [[J00, J01], [J10, J11]]
F = [
    dolfinx.fem.form(F_0, entity_maps=entity_maps),
    dolfinx.fem.form(F_1, entity_maps=entity_maps),
]
b_bc = dolfinx.fem.Function(u_0.function_space)
b_bc.x.array[:] = 0.2
submesh_b.topology.create_connectivity(
    submesh_b.topology.dim - 1, submesh_b.topology.dim
)
bc_b = dolfinx.fem.dirichletbc(
    b_bc, dolfinx.fem.locate_dofs_topological(u_0.function_space, fdim, ft_b.find(2))
)

t_bc = dolfinx.fem.Function(u_1.function_space)
t_bc.x.array[:] = 0.05
submesh_t.topology.create_connectivity(
    submesh_t.topology.dim - 1, submesh_t.topology.dim
)
bc_t = dolfinx.fem.dirichletbc(
    t_bc, dolfinx.fem.locate_dofs_topological(u_1.function_space, fdim, ft_t.find(1))
)
bcs = [bc_b, bc_t]

solver = NewtonSolver(
    F,
    J,
    [u_0, u_1],
    bcs=bcs,
    max_iterations=2,
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
        "pc_factor_mat_solver_type": "mumps",
    },
)
solver.solve(1e-5)

bp = dolfinx.io.VTXWriter(mesh.comm, "u_b.bp", [u_0], engine="BP4")
bp.write(0)
bp.close()
bp = dolfinx.io.VTXWriter(mesh.comm, "u_t.bp", [u_1], engine="BP4")
bp.write(0)
bp.close()

```

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [June 20, 2024, 7:38pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/3 "2024-06-20T19:38:54Z")

</div>

Yes it works thank you so much!

Just a typo

```auto
F_0 = (
    -0.5 * mixed_term((u_b + u_t), v_b, n_b) * dInterface
    - 0.5 * mixed_term(v_b, (u_b / K_b - u_t / K_t), n_b) * dInterface
)

F_1 = (
    +0.5 * mixed_term((u_b + u_t), v_t, n_b) * dInterface
    - 0.5 * mixed_term(v_t, (u_b / K_b - u_t / K_t), n_b) * dInterface
)
F_0 += 2 * gamma / (h_b + h_t) * (u_b / K_b - u_t / K_t) * v_b * dInterface
F_1 += -2 * gamma / (h_b + h_t) * (u_b / K_b - u_t / K_t) * v_t * dInterface

```

---

<div class="post-metadata">

**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, 7:47pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/4 "2024-06-20T19:47:41Z")

</div>

Right, I was a bit too fast 🙂 Corrected in post now:)

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [June 20, 2024, 7:49pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/5 "2024-06-20T19:49:34Z")

</div>

It would be nice to have a demo somewhere to be able to arrive to this by hand

---

<div class="post-metadata">

**Author:** ![RemDelaporteMathurin](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/remdelaportemathurin/32/4543_2.png) [@RemDelaporteMathurin](https://fenicsproject.discourse.group/u/RemDelaporteMathurin)\
**Post date:** [June 20, 2024, 7:51pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/6 "2024-06-20T19:51:01Z")

</div>

![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/a/afb501d3b390d8215127521007a8959630e1f3c0.png)  
For anyone interested this is the produced field

---

<div class="post-metadata">

**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, 8:08pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/7 "2024-06-20T20:08:33Z")

</div>

In the referenced gist, I’ve written out all the different ways of writing up the DG coupling for general DG, which kind of explains the evolution of it.  
I’ll talk to @jpdean about making an official demo out of this (or the linear version of it) when I see him early July:)

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [September 5, 2024, 1:47am UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/8 "2024-09-05T01:47:55Z")

</div>

Have you written a preprint on the variation formulation in this [gist](https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py) or is the prescription of discontinuity by that method supposed to be obvious?

---

<div class="post-metadata">

**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:** [September 5, 2024, 5:57am UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/9 "2024-09-05T05:57:13Z")

</div>

It is a «standard» DG coupling scheme of a Poisson equation (see the dg.py file for a similar DG formulation).

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [September 5, 2024, 7:17pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/10 "2024-09-05T19:17:21Z")

</div>

Thanks for the gist - symmetric interior penalty DG (SIPDG) extra DoFs were becoming a bottleneck for a problem I am solving. Reviewing the gist, I can see the similarities with SIPDG. Do you think the penalty parameter gamma would be affected by mesh quality far away from the interface? I’m trying to ascertain whether I can use bounds on the penalty parameters already proved by [Epshteyn and Riviere](https://www.cmor-faculty.rice.edu/~br1/papers/ERJCAM.pdf). Regardless, I’ll most likely get away with using a “large-enough” penalty term by trial and error if the bounds for this DG are not contained in the bounds by Epshteyn and Riviere.

---

<div class="post-metadata">

**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:** [September 5, 2024, 9:25pm UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/11 "2024-09-05T21:25:05Z")

</div>

I would be suprised if far-away bad cells would affect the DG coupling at the interface itself.  
However, i am aware of cases where boundary effects reduce global convergence (in fluid flow), so there might be similarities.

---

<div class="post-metadata">

**Author:** ![leshinka](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/leshinka/32/5789_2.png) [@leshinka](https://fenicsproject.discourse.group/u/leshinka)\
**Post date:** [September 13, 2024, 12:09am UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/12 "2024-09-13T00:09:35Z")

</div>

I modified the gist by:

1. importing gmsh
2. defining a function that creates a mesh and saves to file “mesh.msh”
3. loading the mesh from file using dolfinx.io.gmshio

When I run the problem, clear the cache `rm -r ~/.cache/fenics` and rerun the problem, I get segmentation error.

```auto
===================================================================================
= BAD TERMINATION OF ONE OF YOUR APPLICATION PROCESSES
= PID 142239 RUNNING AT lmolel-thinkpad-t15-gen-2i
= EXIT CODE: 139
= CLEANING UP REMAINING PROCESSES
= YOU CAN IGNORE THE BELOW CLEANUP MESSAGES
===================================================================================
YOUR APPLICATION TERMINATED WITH THE EXIT STRING: Segmentation fault (signal 11)
This typically refers to a problem with your application.
Please see the FAQ page for debugging suggestions

```

Is this supposed to happen, and how do I go about debugging it?

Here’s the function creating mesh:

```auto
def create_mesh():
    points_left = [
        (0, 0, 0),
        (0, 0.5, 0),
        (0, 1, 0),
    ]

    points_right = [
        (1, 0, 0),
        (1, 0.5, 0),
        (1, 1, 0)
    ]

    gmsh.initialize()
    gmsh.model.add('2d')
    left_points = [gmsh.model.occ.addPoint(*p) for p in points_left]
    right_points = [gmsh.model.occ.addPoint(*p) for p in points_right]

    lines_left = [gmsh.model.occ.addLine(left_points[i], left_points[i+1]) for i in range(2)]
    lines_right = [gmsh.model.occ.addLine(right_points[i], right_points[i+1]) for i in range(2)]
    line_top = gmsh.model.occ.addLine(left_points[-1], right_points[-1])
    line_middle = gmsh.model.occ.addLine(left_points[1], right_points[1])
    line_bottom = gmsh.model.occ.addLine(left_points[0], right_points[0])

    loop_top = gmsh.model.occ.addCurveLoop([lines_left[-1], line_middle, lines_right[-1], line_top])
    top_surf = gmsh.model.occ.addPlaneSurface([loop_top])
    loop_bottom = gmsh.model.occ.addCurveLoop([lines_left[0], line_bottom, lines_right[0], line_middle])
    bottom_surf = gmsh.model.occ.addPlaneSurface([loop_bottom])
    gmsh.model.occ.synchronize()
    gmsh.model.addPhysicalGroup(1, [line_top], 1, "top")
    gmsh.model.addPhysicalGroup(1, [line_bottom], 2, "bottom")
    gmsh.model.addPhysicalGroup(1, [line_middle], 5, "interface")
    gmsh.model.addPhysicalGroup(2, [bottom_surf], 3, "bottom")
    gmsh.model.addPhysicalGroup(2, [top_surf], 4, "top")
    gmsh.model.mesh.generate(2)
    gmsh.write("mesh.msh")
    gmsh.finalize()

```

In the [gist](https://gist.githubusercontent.com/jorgensd/8a5c32f491195e838f5863ca88b27bce/raw/697faabb347481cb1da1f1f7a5b3d2d4ea309e92/problem.py), I commented out lines 196-230 and inserted the following snippet at line 231

```auto
# load mesh from file
comm = MPI.COMM_WORLD
partitioner = dolfinx.mesh.create_cell_partitioner(dolfinx.mesh.GhostMode.shared_facet)
mesh, ct, mt = dolfinx.io.gmshio.read_from_msh("mesh.msh", comm, partitioner=partitioner)
gdim = mesh.geometry.dim
tdim = mesh.topology.dim
fdim = tdim - 1
num_facets_local = (
    mesh.topology.index_map(fdim).size_local + mesh.topology.index_map(fdim).num_ghosts
)

```

I encounter the issue in conda for `fenics-dolfinx==0.8.0.0` and not for local build of development `fenics-dolfinx==0.9.0.0`.

---

<div class="post-metadata">

**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:** [September 13, 2024, 3:36am UTC](https://fenicsproject.discourse.group/t/discontinuity-at-interface-using-mixed-domains/15040/13 "2024-09-13T03:36:26Z")

</div>

Please note that I have addresses several issues Since the 0.8 release when it comes to submeshes,

> <https://github.com/FEniCS/dolfinx/pull/3369>
>
> Currently, collapsing a dofmap with no dofs on the process segfaults due to unsa…fe access of \`dofs\_view.back()\`, as it might be empty.
> 
> This is unlikely to happen to "standard" meshes, but very likely to happen for sub-meshes.
> 
> Thanks to @RemDelaporteMathurin for providing me with a segfaulting code that I could distill down.

> <https://github.com/FEniCS/dolfinx/pull/3361>
>
> Continuation of https://github.com/FEniCS/dolfinx/pull/3260.
> 
> If one has multi…ple disjoint sub-meshes in a variational form, then restricting packing by \`enabled\_coefficients\` per integral type is not sufficient.
> 
> Example follows:
> Divide an interval mesh in three disjoint cell sets (left, center, right).
> If we want to create a coupling from the right submesh to the parent mesh at the interface between center and right),
> we can add an integral \`(u\_parent("+")-u\_parent("-"))\*v\_right("+")\*dS(interface\_rc)\` to our variational form.
> Similarly we could do the same for the left integral: \`(u\_parent("+")-u\_parent("-"))\*v\_left("+")\*dS(interface\_lc)\`.
> 
> However, now both coefficient v\_right and v\_left will be packed for all facets marked with either interface\_rc or interface\_lc, which leads to undefined behavior (segfaults).
> 
> The following PR adds a check after fetching the cell during packing, and only packs coefficients if the cell index is positive.
> 
> Issue reported by @RemDelaporteMathurin when trying to adapt my example: https://gist.github.com/jorgensd/8a5c32f491195e838f5863ca88b27bce#file-problem-py to multiple surfaces.

> <https://github.com/FEniCS/dolfinx/pull/3260>
>
> Pack only entities used in integral type:
> MWE:
> \`\`\`python
> from mpi4py import M…PI
> import dolfinx
> import ufl
> 
> mesh = dolfinx.mesh.create\_unit\_cube(MPI.COMM\_WORLD, 50, 50, 50)
> V = dolfinx.fem.functionspace(mesh, ("Lagrange", 3))
> N = 50
> us = \[dolfinx.fem.Function(V) for \_ in range(N)\]
> for i,u in enumerate(us):
> u.x.array\[:\] = i + 1
> u.name=f"u\_{i}"
> 
> volume\_form = us\[0\] \* ufl.dx
> surface\_form = sum(us\[i\] \* ufl.ds for i in range(N))
> 
> combined\_form = volume\_form + surface\_form
> compiled\_form = dolfinx.fem.form(combined\_form)
> 
> import time
> start = time.perf\_counter()
> coeffs = dolfinx.fem.assemble.pack\_coefficients(compiled\_form.\_cpp\_object)
> end =time.perf\_counter()
> print(f"{dolfinx.common.git\_commit\_hash=} {end-start=:.5e}")
> print(coeffs\[(dolfinx.fem.IntegralType.cell, -1)\])
> \`\`\`
> On the main branch, all coefficients are packed for the cell integral, even if only the first coefficient is used in that integral.
> Runtime on this branch:
> \`\`\`bash
> dolfinx.common.git\_commit\_hash='002e1fa9360995f9c71414a683e7ce618d6101e9' end-start=1.70250e+00
> \`\`\`
> on main
> \`\`\`bash
> dolfinx.common.git\_commit\_hash='64e35310085683f4f61804e70b8b58b9bb8653ae' end-start=2.89335e+00
> \`\`\`
> It would also resolve #3256 and simplify some logic in: #3224