Using 2 disconnected meshes and connecting them with jump condition

General
9

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

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 and post2, 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.