TypeError in 3D electromagnetics problem

General
1

Hi,

I am trying to solve a 3D electromagnetic wave propagation problem in which the incident wave is excited from a port boundary at the top (z=fov_z/2) and the same is absorbed by another port boundary at the bottom (z=-fov_z/2). I cannot figure my way on the following:

  1. How to define nxE=0 boundary condition where I want? I mean, how do we code that only the dofs corresponding to the tangential components of the field to be zero and not the normal?
  2. An error message I cannot crack when I run the code.

Thanks in advance for the help. Below is my MWE which I am running on dolfinx 0.6.0:

import dolfinx, mpi4py, gmsh, ufl, petsc4py
from dolfinx.io import gmshio
import numpy as np
import typing

#%% Input parameters
fov_x, fov_y, fov_z = 1, 1, 2 # in um
t_film = 0.1 # film thickness
dpml = 0.1 # thickness of pmls

port_bnd = (6, 5) # incident and transmitted port

wlen = 1 # in um
U0 = 1 # incident amplitude
theta = 0

#%%

#%% Create geometry

def create_geom():
    gmsh.initialize()
    gmsh.option.setNumber("General.Terminal", 0) # disable output messages
    gmsh.clear()
    gdim = 3
    model_rank = 0

    box1 = gmsh.model.occ.addBox(-fov_x/2, -fov_y/2, -fov_z/2, fov_x, fov_y, fov_z)
    gmsh.model.occ.synchronize()
    
    # number all domains
    all_doms = gmsh.model.getEntities(gdim)
    for j, edge in enumerate(all_doms):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1)  # create the main group/node

    # number all boundaries
    all_edges = gmsh.model.getEntities(gdim - 1)
    for j, edge in enumerate(all_edges):
        gmsh.model.addPhysicalGroup(edge[0], [edge[1]], j + 1)  # create the main group/node

    gmsh.model.mesh.generate(gdim)
    gmsh.model.mesh.refine()
    # gmsh.write('geom/plane_wave_mwe.msh')
    model_rank = 0
    mesh, ct, ft = gmshio.model_to_mesh(gmsh.model, mpi4py.MPI.COMM_WORLD, model_rank, gdim)
    gmsh.fltk.run() # visible inspection, for debugging
    gmsh.finalize()
    return mesh, ct, ft

mesh, ct, ft = create_geom() # create fenics geometry, cell tags and facet tags
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

#%%

elem = ufl.FiniteElement('Nedelec 1st kind H(curl)', mesh.ufl_cell(), degree=2)
Omega = dolfinx.fem.FunctionSpace(mesh, elem)

k0 = dolfinx.fem.Constant(mesh, 2*np.pi/wlen) # vacuum wave number
# #%%

U = ufl.TrialFunction(Omega)
Ut = ufl.TestFunction(Omega)
dx = ufl.Measure("dx", mesh, subdomain_data=ct) # for selective integration in domains
p_param = 1 # to be modified later
q_param = 1 # to be modified later

# weak form in the volume
F = (p_param*ufl.inner(ufl.curl(U), ufl.curl(Ut)) - k0**2*q_param*ufl.inner(U, Ut))*ufl.dx

#%% boundary conditions
class PlaneWave:
    def __init__(self, U0:tuple, k0:float, n:complex=1, theta:float=0):
        self.U0 = U0 # amplitude
        self.theta = theta  # incident angle
        self.k0 = k0  # vacuum wavevector
        self.n = n  # refractive index
        
    def eval(self, x: np.typing.NDArray[np.float64]) -> typing.Tuple[
            np.typing.NDArray[np.complex128],
            np.typing.NDArray[np.complex128]]:

        kx = self.k0 * self.n * np.sin(self.theta)
        ky = 0
        kz = self.k0 * self.n * np.cos(self.theta)
        
        phi = kx * x[0] + ky * x[1] + kz * x[2]

        return (self.U0[0] * np.exp(1j * phi), self.U0[1] * np.exp(1j * phi), self.U0[2] * np.exp(1j * phi))

n = ufl.FacetNormal(mesh) # boundary normal vector
ds = ufl.Measure("ds", domain=mesh, subdomain_data=ft) # boundary differential element

kxi = dolfinx.fem.Constant(mesh, 0.0)
kzi = dolfinx.fem.Constant(mesh, np.sqrt(k0.value**2-kxi.value**2)) # incident medium
kzt = dolfinx.fem.Constant(mesh, np.sqrt(k0.value**2-kxi.value**2)) # transmitted medium

ki = ufl.as_vector((kxi, 0, -kzi)) # incident, going in -z direction
kr = ufl.as_vector((kxi, 0, kzi)) # reflected, going in +z direction
kt = ufl.as_vector((kxi, 0, -kzt)) # transmitted, going in -z direction

# U_inc = ufl.as_vector((0, U0, 0))
pw_inc = PlaneWave((0,U0,0), k0.value, 1, 0)
U_inc = dolfinx.fem.Function(Omega)
U_inc.interpolate(pw_inc.eval)

bc_inc_abs = ufl.inner(-1j*kzi*U, Ut)*ds(port_bnd[0]) # incident port
bc_tr_abs = ufl.inner(-1j*kzt*U, Ut)*ds(port_bnd[1]) # transmission port
bc_inc_exc = ufl.inner(2j*kzi*U_inc, Ut)*ds(port_bnd[0]) # incident port excitation

F += bc_inc_abs
F += bc_tr_abs
F += bc_inc_exc

a = ufl.lhs(F)
L = ufl.rhs(F)

# dirichlet: set the polarization componets x,z zero
# for p-pol, it amounts to choosing PMC while for s-polarization PEC
# QUESTION: Is this correct way of doing nxE = 0?
def perfect_cond_bnd(x:np.typing.NDArray[np.float64]):
    '''identify pec/pmc planes for s/p polarization'''
    return np.isclose(np.abs(x[1]), fov_y/2)
gdim=3
perfect_cond_facets = dolfinx.mesh.locate_entities_boundary(mesh, gdim-1, perfect_cond_bnd)
perfect_cond_dofs = dolfinx.fem.locate_dofs_topological(Omega, mesh.topology.dim-1, perfect_cond_facets)
ubc = dolfinx.fem.Function(Omega)
ubc.x.set(0+0j)
bc_perfect_cond = dolfinx.fem.dirichletbc(ubc, perfect_cond_dofs)

#%% solve
Lin_system = dolfinx.fem.petsc.LinearProblem(a, L, bcs = [bc_perfect_cond], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
U_sol = Lin_system.solve()