Hi,

I am solving vector Helmholtz equation in 3D and need perfect electric conductors (PEC) at both external and internal boundaries in my geometry. For external boundaries, I am able to realize it with Nitche’s method. However, I am unsure about internal boundaries due to field discontinuities. For now, I just add Nitche’s term on both sides of the edges to the variational form.

Nonetheless, I check my approach with a simple toy model. I divide a vacuum box into two halves separated by a PEC and launch a plane wave in only one. If the internal PEC works reliably, half of the volume should be dark at noise levels. But I get observable amplitudes in the supposedly dark box too. Below is the paraview output:

I will appreciate if knowledgeable users could help me understand how the Nitche’s approach is modified for internal facets or point me to literature where I might find a solution. I have the book of Jian-Ming Jin but it seems silent on this.

Below is the simple MWE which produced the result shown above. To avoid mesh tags from changing (hopefully), the gmsh mesh file can be downloaded from here.

```
import dolfinx, petsc4py, mpi4py, gmsh, ufl
from dolfinx.io import gmshio
import numpy as np
from dolfinx.fem.petsc import LinearProblem
import time
#%% Load mesh
gmsh.initialize()
gmsh.option.setNumber("General.Terminal", 0) # disable output messages
gmsh.clear()
gmsh.merge('mwe.msh')
model_rank = 0
gdim = 3
mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
#%% input paramters
port_bnd_idx = (11, 10) # incident, transmitted
pec_int_bnd_idx = (2, )
pec_ext_bnd_idx = (9, 4)
k0 = 2*np.pi/0.5
fe_order = 2
#%% Variational form
Omega = dolfinx.fem.functionspace(mesh, ("N1curl", fe_order))
U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)
F = (ufl.inner(ufl.curl(U), ufl.curl(Ut)) - k0**2*ufl.inner(U, Ut))*ufl.dx
F += -ufl.inner(dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType((0.0, 0.0,0.0))), Ut)*ufl.dx # source term
#%% Port boundary at the top (excitation) and bottom surfaces
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft) # boundary differential element
n_3d = ufl.FacetNormal(mesh)
ki_3d = ufl.as_vector((0, 0, -k0)) # incident wave
kr_3d = ufl.as_vector((0, 0, k0)) # reflected wave, if any
kt_3d = ufl.as_vector((0, 0, -k0)) # transmitted wave
# port boundary at the incident port
# reflected wave
F_bnd_term = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*kr_3d, U)), Ut)*ds(port_bnd_idx[0])
F += F_bnd_term
# excitation
E_inc = ufl.as_vector((0, 1, 0)) # y-polarized
F_bnd_term_exc = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*(ki_3d-kr_3d), E_inc)), Ut)*ds(port_bnd_idx[0])
F += F_bnd_term_exc
# port boundary at the transmitted port
F_bc_tr = ufl.inner(ufl.cross(n_3d, ufl.cross(1j*kt_3d, U)), Ut)*ds(port_bnd_idx[1])
F += F_bc_tr
#%% PEC boundaries using Nitche's method
# hE = 2 * ufl.Circumradius(mesh)
hE = ufl.CellDiameter(mesh)
alpha = dolfinx.fem.Constant(mesh, petsc4py.PETSc.ScalarType(10))
U_pec = ufl.as_vector((U[0],0,U[2])) # Non-zero components will be forced zero
F_pec_ext = ufl.inner(ufl.cross(n_3d, ufl.curl(U)), Ut)*ds(pec_ext_bnd_idx) + alpha/hE*ufl.inner(U_pec, Ut)*ds(pec_ext_bnd_idx) - \
ufl.inner(U_pec, ufl.cross(n_3d, ufl.curl(Ut)))*ds(pec_ext_bnd_idx)
dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft) # boundary differential element
F_pec_int_pos = ufl.inner(ufl.cross(n_3d, ufl.curl(U))('+'), Ut('+'))*dS(pec_int_bnd_idx) + alpha/hE('+')*ufl.inner(U_pec('+'), Ut('+'))*dS(pec_int_bnd_idx) - \
ufl.inner(U_pec('+'), ufl.cross(n_3d, ufl.curl(Ut))('+'))*dS(pec_int_bnd_idx)
F_pec_int_neg = ufl.inner(ufl.cross(n_3d, ufl.curl(U))('-'), Ut('-'))*dS(pec_int_bnd_idx) + alpha/hE('-')*ufl.inner(U_pec('-'), Ut('-'))*dS(pec_int_bnd_idx) - \
ufl.inner(U_pec('-'), ufl.cross(n_3d, ufl.curl(Ut))('-'))*dS(pec_int_bnd_idx)
F += F_pec_ext
F += F_pec_int_pos
F += F_pec_int_neg
a = ufl.lhs(F)
L = ufl.rhs(F)
#%% Solve
U_sol = dolfinx.fem.Function(Omega)
petsc_options = {"ksp_type": "preonly", "pc_type": "lu", "pc_factor_mat_solver_type":"mumps"}
Lin_system = LinearProblem(a, L, bcs = [], u=U_sol, petsc_options=petsc_options)
solver = Lin_system.solver
t1 = time.perf_counter()
Lin_system.solve()
t2 = time.perf_counter()
print("Solved in {}s".format(t2-t1))
Omega_sol = dolfinx.fem.functionspace(mesh, ("DG", fe_order, (mesh.topology.dim,)))
U_sol_dg = dolfinx.fem.Function(Omega_sol)
U_sol_dg.interpolate(U_sol)
with dolfinx.io.VTXWriter(mpi4py.MPI.COMM_WORLD, "E_field.bp", [U_sol_dg], engine="BP4") as vtx:
vtx.write(0.0)
```