# Using 2 disconnected meshes and connecting them with jump condition

**URL:** <https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304>\
**Category:** General\
**Created:** [January 31, 2023, 12:41pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304 "2023-01-31T12:41:49Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [January 31, 2023, 12:41pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/1 "2023-01-31T12:41:49Z")

</div>

Hi there,

I want to impose a jump condition between two separate meshes like below;

 ![Screenshot from 2023-01-31 12-36-45](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/8/8d93eecea95a0e79c35a309ea7265bc6fe2d69d0.png)

Is it possible in dolfinx? Should define two different meshes or mixed elements for one mesh?

I did try using 2 seperate mesh files and when I read subdomains and facets\_tags, I get;

```auto
File "main.py", line 33, in <module>
    matrices = PassiveFlame2(mesh1, mesh2, boundary_conditions, c, degree=degree)
  File "/home/ee331/Dev/Acoustics-Dev/LOTAN/HelmholtzX/MP/Javier/passive_flame_x.py", line 280, in __init__
    self.a_form = form(-self.c**2 * inner(grad(self.u1), grad(self.v1))*self.dx1(1)-self.c**2 * inner(grad(self.u2), grad(self.v2))*self.dx2(2))
  File "/home/ee331/Dev/Venvs/v042c/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 139, in form
    return _create_form(form)
  File "/home/ee331/Dev/Venvs/v042c/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 134, in _create_form
    return _form(form)
  File "/home/ee331/Dev/Venvs/v042c/lib/python3.8/dist-packages/dolfinx/fem/forms.py", line 102, in _form
    subdomains, = list(sd.values()) # Assuming single domain
ValueError: too many values to unpack (expected 1)

```

I think dolfinx can only handle single domain. But wonder is there any other way of implementing this?

Thank you for your answers in advance.

---

<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:** [February 2, 2023, 1:48pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/2 "2023-02-02T13:48:08Z")

</div>

There is no way with `ufl` to prescribe such a variational form.

Depending on what the `jump` condition between the meshes is supposed to do, you could use `dolfinx_mpc`.

---

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [February 2, 2023, 2:12pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/3 "2023-02-02T14:12:53Z")

</div>

Can I still use separate meshes, if I map these 2 boundaries using `dolfinx_mpc`?

---

<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:** [February 2, 2023, 2:14pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/4 "2023-02-02T14:14:11Z")

</div>

See for instance: [dolfinx\_mpc/demo\_elasticity\_disconnect\_2D.py at main · jorgensd/dolfinx\_mpc · GitHub](https://github.com/jorgensd/dolfinx_mpc/blob/main/python/demos/demo_elasticity_disconnect_2D.py)  
Depending on your `jump`-condition, you can do it in `dolfinx_mpc`.

---

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [February 2, 2023, 3:44pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/5 "2023-02-02T15:44:14Z")

</div>

I am solving an eigenvalue problem. Weak form of my equation is;

-\int\_\Omega s^2 \nabla u \cdot \nabla \hat{v} d\Omega + \int\_{\partial \Omega} s^2 (\nabla \hat{u} \cdot \mathbf{n}) v dS + \omega^2 \int\_\Omega u \cdot \hat{v} d\Omega = 0

where

\nabla \hat{u} \cdot \mathbf{n} = [\hat{u}\_2-\hat{u}\_1] for the mesh below;

 ![Screenshot from 2023-02-02 15-42-28](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6792a03efb1e0b4eaea0bf3e83e06e6830af4570.png)

I couldn’t understand how I impose the pressure difference on boundaries. Here is my MWE;

```auto
import gmsh
import numpy as np
from typing import List, Tuple
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace,
                         VectorFunctionSpace, locate_dofs_geometrical,
                         locate_dofs_topological,form)
from dolfinx.io import XDMFFile
from dolfinx.mesh import locate_entities_boundary
from dolfinx_mpc import LinearProblem, MultiPointConstraint, assemble_matrix
from dolfinx_mpc.utils import (create_point_to_point_constraint,
                               gmsh_model_to_mesh, rigid_motions_nullspace)
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (Identity, Measure, SpatialCoordinate, TestFunction,
                 TrialFunction, grad, inner, sym, tr)
from slepc4py import SLEPc

# Mesh parameters for creating a mesh consisting of two disjoint rectangles
right_tag = 1
left_tag = 2

gmsh.initialize()
if MPI.COMM_WORLD.rank == 0:
    gmsh.clear()

    left_rectangle = gmsh.model.occ.addRectangle(0, 0, 0, 1, 1)
    right_rectangle = gmsh.model.occ.addRectangle(1.1, 0, 0, 1, 1)

    gmsh.model.occ.synchronize()

    # Add physical tags for volumes
    gmsh.model.addPhysicalGroup(1, [8], tag=2)
    gmsh.model.setPhysicalName(1, 2, "Downstream")
    gmsh.model.addPhysicalGroup(1, [2], tag=1)
    gmsh.model.setPhysicalName(1, 1, "Upstream")

    gmsh.model.addPhysicalGroup(2, [right_rectangle], tag=right_tag)
    gmsh.model.setPhysicalName(2, right_tag, "Right square")
    gmsh.model.addPhysicalGroup(2, [left_rectangle], tag=left_tag)
    gmsh.model.setPhysicalName(2, left_tag, "Left square")
    gmsh.option.setNumber("Mesh.CharacteristicLengthMin", 1)
    gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 1)
    gmsh.option.setNumber("Mesh.MeshSizeMax", 0.1)
    # Generate mesh
    gmsh.model.mesh.generate(2)

    # gmsh.option.setNumber("General.Terminal", 1)
    gmsh.model.mesh.optimize("Netgen")
    # gmsh.fltk.run()

mesh, ct,ft = gmsh_model_to_mesh(gmsh.model, cell_data=True, facet_data=True, gdim=2)
gmsh.clear()
gmsh.finalize()
MPI.COMM_WORLD.barrier()

with XDMFFile(mesh.comm, "test.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh)

V = FunctionSpace(mesh, ("Lagrange", 1))
tdim = mesh.topology.dim
fdim = tdim - 1

DG0 = FunctionSpace(mesh, ("DG", 0))
left_cells = ct.indices[ct.values == left_tag]
right_cells = ct.indices[ct.values == right_tag]
left_dofs = locate_dofs_topological(DG0, tdim, left_cells)
right_dofs = locate_dofs_topological(DG0, tdim, right_cells)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
dx = Measure("dx", domain=mesh, subdomain_data=ct)
ds = Measure("ds", domain=mesh, subdomain_data=ft)

def gather_dof_coordinates(V, dofs):
    """
    Distributes the dof coordinates of this subset of dofs to all processors
    """
    x = V.tabulate_dof_coordinates()
    local_dofs = dofs[dofs < V.dofmap.index_map.size_local * V.dofmap.index_map_bs]
    coords = x[local_dofs]
    num_nodes = len(coords)
    glob_num_nodes = MPI.COMM_WORLD.allreduce(num_nodes, op=MPI.SUM)
    recvbuf = None
    if MPI.COMM_WORLD.rank == 0:
        recvbuf = np.zeros(3 * glob_num_nodes, dtype=np.float64)
    sendbuf = coords.reshape(-1)
    sendcounts = np.array(MPI.COMM_WORLD.gather(len(sendbuf), 0))
    MPI.COMM_WORLD.Gatherv(sendbuf, (recvbuf, sendcounts), root=0)
    glob_coords = MPI.COMM_WORLD.bcast(recvbuf, root=0).reshape((-1, 3))
    return glob_coords

facets_r = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 1))
dofs_r = locate_dofs_topological(V, fdim, facets_r)
facets_l = locate_entities_boundary(mesh, fdim, lambda x: np.isclose(x[0], 1.1))
dofs_l = locate_dofs_topological(V, fdim, facets_l)

# Given the local coordinates of the dofs, distribute them on all processors
nodes = [gather_dof_coordinates(V, dofs_r), gather_dof_coordinates(V, dofs_l)]
pairs = []
for x in nodes[0]:
    for y in nodes[1]:
        if np.isclose(x[1], y[1]):
            pairs.append([x, y])
            break

mpc = MultiPointConstraint(V)
for i, pair in enumerate(pairs):
    sl, ms, co, ow, off = create_point_to_point_constraint(
        V, pair[0], pair[1])
    mpc.add_constraint(V, sl, ms, co, ow, off)
mpc.finalize()

s = 343
a_form = -s**2 * inner(grad(u), grad(v))*dx

e_form = s**2 * (inner(u, v)*ds(2)-inner(u, v)*ds(1)) # THIS MIGHT BE WRONG

mass_form = form(a_form+e_form)

A = assemble_matrix(mass_form, mpc)
A.assemble()

stiffnes_form = form(inner(u , v) * dx)

C = assemble_matrix(stiffnes_form, mpc)
C.assemble()

def print0(string: str):
    """Print on rank 0 only"""
    if MPI.COMM_WORLD.rank == 0:
        print(string)

def solve_GEP_shiftinvert(A: PETSc.Mat, B: PETSc.Mat,
                          problem_type: SLEPc.EPS.ProblemType = SLEPc.EPS.ProblemType.GNHEP,
                          solver: SLEPc.EPS.Type = SLEPc.EPS.Type.KRYLOVSCHUR,
                          nev: int = 10, tol: float = 1e-7, max_it: int = 10,
                          target: float = 0.0, shift: float = 0.0) -> SLEPc.EPS:

    # Build an Eigenvalue Problem Solver object
    EPS = SLEPc.EPS()
    EPS.create(comm=MPI.COMM_WORLD)
    A=-A
    EPS.setOperators(A, B)
    EPS.setProblemType(problem_type)
    # set the number of eigenvalues requested
    EPS.setDimensions(nev=nev)
    # Set solver
    EPS.setType(solver)
    # set eigenvalues of interest
    EPS.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE)
    EPS.setTarget(target) # sorting
    # set tolerance and max iterations
    EPS.setTolerances(tol=tol, max_it=max_it)
    # Set up shift-and-invert
    # Only work if 'whichEigenpairs' is 'TARGET_XX'
    ST = EPS.getST()
    ST.setType(SLEPc.ST.Type.SINVERT)
    ST.setShift(shift)
    EPS.setST(ST)
    # set monitor
    it_skip = 1
    EPS.setMonitor(lambda eps, it, nconv, eig, err:
                   monitor_EPS_short(eps, it, nconv, eig, err, it_skip))
    # parse command line options
    EPS.setFromOptions()
    # Display all options (including those of ST object)
    # EPS.view()
    EPS.solve()
    EPS_print_results(EPS)
    return EPS

def EPS_get_spectrum(EPS: SLEPc.EPS,
                     mpc: MultiPointConstraint) -> Tuple[List[complex], List[PETSc.Vec], List[PETSc.Vec]]:
    """ Retrieve eigenvalues and eigenfunctions from SLEPc EPS object.
        Parameters
        ----------
        EPS
           The SLEPc solver
        mpc
           The multipoint constraint

        Returns
        -------
            Tuple consisting of: List of complex converted eigenvalues,
             lists of converted eigenvectors (real part) and (imaginary part)
        """
    # Get results in lists
    eigval = [EPS.getEigenvalue(i) for i in range(EPS.getConverged())]
    eigvec_r = list()
    eigvec_i = list()
    V = mpc.function_space
    for i in range(EPS.getConverged()):
        vr = Function(V)
        vi = Function(V)

        EPS.getEigenvector(i, vr.vector, vi.vector)
        eigvec_r.append(vr)
        eigvec_i.append(vi) # Sort by increasing real parts
    idx = np.argsort(np.real(np.array(eigval)), axis=0)
    eigval = [eigval[i] for i in idx]
    eigvec_r = [eigvec_r[i] for i in idx]
    eigvec_i = [eigvec_i[i] for i in idx]
    return (eigval, eigvec_r, eigvec_i)

def monitor_EPS_short(EPS: SLEPc.EPS, it: int, nconv: int, eig: list, err: list, it_skip: int):
    """
    Concise monitor for EPS.solve().

    Parameters
    ----------
    eps
        Eigenvalue Problem Solver class.
    it
       Current iteration number.
    nconv
       Number of converged eigenvalue.
    eig
       Eigenvalues
    err
       Computed errors.
    it_skip
        Iteration skip.

    """
    if (it == 1):
        print0(' ******************************')
        print0(' ***SLEPc Iterations...***')
        print0(' ******************************')
        print0("Iter. | Conv. | Max. error")
        print0(f"{it:5d} | {nconv:5d} | {max(err):1.1e}")
    elif not it % it_skip:
        print0(f"{it:5d} | {nconv:5d} | {max(err):1.1e}")

def EPS_print_results(EPS: SLEPc.EPS):
    """ Print summary of solution results. """
    print0("\n ******************************")
    print0(" ***SLEPc Solution Results***")
    print0(" ******************************")
    its = EPS.getIterationNumber()
    print0(f"Iteration number: {its}")
    nconv = EPS.getConverged()
    print0(f"Converged eigenpairs: {nconv}")

    if nconv > 0:
        # Create the results vectors
        vr, vi = EPS.getOperators()[0].createVecs()
        print0("\nConverged eigval. Error ")
        print0("----------------- -------")
        for i in range(nconv):
            k = EPS.getEigenpair(i, vr, vi)
            error = EPS.computeError(i)
            if k.imag != 0.0:
                print0(f" {k.real:2.2e} + {k.imag:2.2e}j {error:1.1e}")
            else:
                pad = " " * 11
                print0(f" {k.real:2.2e} {pad} {error:1.1e}")

target_eig = (170*2*np.pi)**2
EPS = solve_GEP_shiftinvert(A, C, problem_type=SLEPc.EPS.ProblemType.GHEP,
                            solver=SLEPc.EPS.Type.KRYLOVSCHUR,
                            nev=10, tol=1e-7, max_it=10,
                            target=target_eig, shift=1.5)
(eigval, eigvec_r, eigvec_i) = EPS_get_spectrum(EPS, mpc)

for i in range(len(eigval)):
    mpc.backsubstitution(eigvec_r[i].vector)
    mpc.backsubstitution(eigvec_i[i].vector)

print0(f"Computed eigenvalues:\n {np.around(eigval,decimals=2)}")
print0(f"Computed eigenfrequencies:\n {np.sqrt(eigval)/(2*np.pi)}")

# Save first eigenvector
with XDMFFile(MPI.COMM_WORLD, "results/eigenvector.xdmf", "w", encoding=XDMFFile.Encoding.HDF5) as xdmf:
    xdmf.write_mesh(mesh)
    xdmf.write_function(eigvec_r[12])

```

I get this results now;

```auto
Computed eigenvalues:
Computed eigenvalues:
 [-1.00000000e+00+0.j -1.00000000e+00+0.j -1.00000000e+00+0.j
 -1.00000000e+00+0.j -1.00000000e+00+0.j -1.00000000e+00+0.j
 -1.00000000e+00+0.j -1.00000000e+00+0.j -1.00000000e+00+0.j
 -1.00000000e+00+0.j -1.00000000e+00+0.j -0.00000000e+00+0.j
  2.90732620e+05+0.j 1.16776004e+06+0.j 1.16806494e+06+0.j
  1.46205284e+06+0.j]
Computed eigenfrequencies:
 [ 0. +1.59154943e-01j 0. +1.59154943e-01j
   0. +1.59154943e-01j 0. +1.59154943e-01j
   0. +1.59154943e-01j 0. +1.59154943e-01j
   0. +1.59154943e-01j 0. +1.59154943e-01j
   0. +1.59154943e-01j 0. +1.59154943e-01j
   0. +1.59154943e-01j 0. +7.73230399e-07j
  85.81575235+0.00000000e+00j 171.98752123+0.00000000e+00j
 172.009972 +0.00000000e+00j 192.44279967+0.00000000e+00j]

```

These values makes sense, but I wonder whether my implementation is correct or not?

---

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [February 3, 2023, 11:14am UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/6 "2023-02-03T11:14:42Z")

</div>

I cannot get damping (positive imaginary part for eigenvalues) for the eigenvalues. Originally, my damping term is;

\int\_{\partial \Omega} s^2 (\nabla \hat{u} \cdot \mathbf{n} v dS) = \int\_{\partial \Omega} s^2 \frac{K\_R}{d^2}[\hat{u}\_2-\hat{u}\_1]v dS

where

K\_R = 2a\frac{(\pi/2 I\_1(St) e^{-St})-1jK\_1(St)sinh(Sh)}{St((\pi/2 I\_1(St) e^{-St})+1jK\_1(St)cosh(Sh))}, I\_1 and K\_1 are Bessel functions and St = wa/U.

My implementation is;

```auto
c = 347

a = 3*1E-3
d = 35*1E-3
U = 5
omega = (382.9)*2*np.pi
St = omega*a/U
num = (np.pi/2)*jv(1,St)*np.exp(-St)-1j*yv(1,St)*np.sinh(St)
denom = St*((np.pi/2)*jv(1,St)*np.exp(-St)+1j*yv(1,St)*np.cosh(St)) 
expression = 1+num/denom
K_r = 2*a*expression

e_form = -c **2 * K_r/d** 2 * (inner(u, v)*ds(1)-inner(u, v)*ds(2))

```

Is it something mpc cannot handle @dokken ? I am using v0.4.1 of dolfinx\_mpc and v0.4.2 of dolfinx.

---

<div class="post-metadata">

**Author:** ![IngeborgGjerde](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ingeborggjerde/32/4504_2.png) [@IngeborgGjerde](https://fenicsproject.discourse.group/u/IngeborgGjerde)\
**Post date:** [February 4, 2023, 2:23pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/7 "2023-02-04T14:23:54Z")

</div>

Alternatively, you could try adapting the code from [Mirok/emi-book-fem](https://github.com/MiroK/emi-book-fem).  
In particular, take a look at [this demo](https://github.com/MiroK/emi-book-fem/blob/master/emi/fem/prime_multi.py) that shows how to implement Poisson equations on two domains with some coupling condition between them.

The model problem and equations this repo solves are introduced in Chapter 5 of the EMI book: [Solving the EMI Equations using Finite Element Methods | SpringerLink](https://link.springer.com/chapter/10.1007/978-3-030-61157-6_5)

---

<div class="post-metadata">

**Author:** ![Massimiliano\_Leoni](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/massimiliano_leoni/32/584_2.png) [@Massimiliano\_Leoni](https://fenicsproject.discourse.group/u/Massimiliano_Leoni)\
**Post date:** [February 6, 2023, 9:55am UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/8 "2023-02-06T09:55:51Z")

</div>

Implementing your jump condition with dolfinx\_mpc would be quite tricky, I feel.  
You might be better off with a fixed-point iteration, similar to a domain-decomposition method with Robin-type coupling, in which you solve on one domain at a time until you converge.

---

<div class="post-metadata">

**Author:** ![Ekrem\_Ekici](https://avatars.discourse-cdn.com/v4/letter/e/977dab/32.png) [@Ekrem\_Ekici](https://fenicsproject.discourse.group/u/Ekrem_Ekici)\
**Post date:** [February 17, 2023, 4:43pm UTC](https://fenicsproject.discourse.group/t/using-2-disconnected-meshes-and-connecting-them-with-jump-condition/10304/9 "2023-02-17T16:43:15Z")

</div>

I think `dolfinx_mpc` is not suitable for this problem. DG elements could be more useful. I have implemented another simple MWE;

```auto
import gmsh
from mpi4py import MPI
from dolfinx.io import gmshio
from dolfinx.fem import FunctionSpace, form
from dolfinx.fem.petsc import assemble_matrix
from ufl import TestFunction,TrialFunction,Measure,inner,grad,jump
from scipy.special import jv,yv
import numpy as np
from slepc4py import SLEPc

gmsh.initialize()
gdim = 2
mesh_comm = MPI.COMM_WORLD
model_rank = 0
if model_rank == 0:
    gmsh.clear()
    gmsh.option.setNumber("General.Terminal", 0)
    r_inner, r_outer = 0.2 , 0.25
    c1 = gmsh.model.occ.addCircle(0, 0, 0, r_inner)
    cl1 = gmsh.model.occ.addCurveLoop([c1], 1)
    s1 = gmsh.model.occ.addPlaneSurface([cl1])
    c2 = gmsh.model.occ.addCircle(0, 0, 0, r_outer)
    cl2= gmsh.model.occ.addCurveLoop([c2], 2)
    s2 = gmsh.model.occ.addPlaneSurface([cl2], 2)
    gmsh.model.occ.cut([(2,s2)], [(2,s1)],3,True,False)
    gmsh.model.occ.synchronize()
    # Boundary Tags
    gmsh.model.addPhysicalGroup(1, [1], 1)
    gmsh.model.addPhysicalGroup(1, [2], 2)
    # Surface Tags
    gmsh.model.addPhysicalGroup(2, [1], tag=1) # Inner Circle Boundary
    gmsh.model.addPhysicalGroup(2, [3], tag=2) # Outer Circle Boundary
    lc = 0.02 # 0.01
    gmsh.model.occ.synchronize()
    gmsh.option.setNumber("Mesh.MeshSizeMax", lc)
    gmsh.model.mesh.generate(2)
    gmsh.model.mesh.optimize("Netgen")
    # gmsh.fltk.run()

mesh, ct,ft = gmshio.model_to_mesh(gmsh.model, mesh_comm, model_rank, gdim=gdim)
gmsh.clear()
gmsh.finalize()
MPI.COMM_WORLD.barrier()

V = FunctionSpace(mesh, ("DG", 1))
u = TrialFunction(V)
v = TestFunction(V)
dx = Measure('dx', domain=mesh, subdomain_data=ct)
dS = Measure('dS', domain=mesh, subdomain_data=ft)

c = 347
form1 = -c**2 * inner(grad(u), grad(v))*dx

a = 3*1E-3
d = 35*1E-3
U = 5
omega = (382.9)*2*np.pi
St = omega*a/U
num = (np.pi/2)*jv(1,St)*np.exp(-St)-1j*yv(1,St)*np.sinh(St)
denom = St*((np.pi/2)*jv(1,St)*np.exp(-St)+1j*yv(1,St)*np.cosh(St)) 
expression = 1+num/denom
K_r = 2*a*expression
form2 = -c **2 * K_r/d** 2 * inner(jump(u), jump(v))*dS(1)

mass_form = form(form1+form2)
A = assemble_matrix(mass_form)
A.assemble()

stiffnes_form = form(inner(u , v) * dx)
C = assemble_matrix(stiffnes_form)
C.assemble()

def results(E):
    if MPI.COMM_WORLD.Get_rank()==0:
        print()
        print(" ******************************")
        print(" ***SLEPc Solution Results***")
        print(" ******************************")
        print()
        its = E.getIterationNumber()
        print("Number of iterations of the method: %d" % its)
        eps_type = E.getType()
        print("Solution method: %s" % eps_type)
        nev, ncv, mpd = E.getDimensions()
        print("Number of requested eigenvalues: %d" % nev)
        tol, maxit = E.getTolerances()
        print("Stopping condition: tol=%.4g, maxit=%d" % (tol, maxit))
        nconv = E.getConverged()
        print("Number of converged eigenpairs %d" % nconv)
        A = E.getOperators()[0]
        vr, vi = A.createVecs()
        if nconv > 0:
            print()
        for i in range(nconv):
            k = E.getEigenpair(i, vr, vi)
            print("%15f, %15f" % (k.real, k.imag))
        print()

def eps_solver(A, C, target, nev, two_sided=False, print_results=False):
    E = SLEPc.EPS().create(MPI.COMM_WORLD)
    C = - C
    E.setOperators(A, C)
    # spectral transformation
    st = E.getST()
    st.setType('sinvert')
    # E.setKrylovSchurPartitions(1) # MPI.COMM_WORLD.Get_size()
    E.setTarget(target)
    E.setWhichEigenpairs(SLEPc.EPS.Which.TARGET_MAGNITUDE) # TARGET_REAL or TARGET_IMAGINARY
    E.setTwoSided(two_sided)
    E.setDimensions(nev, SLEPc.DECIDE)
    E.setTolerances(1e-15)
    E.setFromOptions()
    E.solve()
    if print_results and MPI.COMM_WORLD.rank == 0:
        results(E)
    return E

target = 400*2*np.pi
nev=10
E = eps_solver(A, C, target**2, nev=nev, print_results= True)
for i in range(nev):
    eig = E.getEigenvalue(i)
    omega = np.sqrt(eig)
    print(omega/(2*np.pi))

```

I should get the eigenvalues around 380Hz with a negative complex part (around -20j). However, I still get **0** solution vectors and eigenvalues, if I use DG elements. According to [post1](https://fenicsproject.discourse.group/t/prescribe-discontinuity-on-internal-boundary/3365/4) and [post2](https://fenicsproject.discourse.group/t/imposing-a-discontinuity-in-the-solution-on-an-internal-interface/10112/7), I might need to manipulate my weak form due to the DG elements usage. But I am not sure how to do it. Is there anybody who can help?

Thank you for your responses in advance.