Solving for a 3D vector field over a 2D mesh with N1curl elements

variational formulation
1

Hi all,

I would like to use fenics for solving this electromagnetic problem: finding the scattered field due to an electromagnetic plane wave (linearly polarized along an arbitrary axis) that hits a metallic cylinder, whose axis of symmetry is perpendicular to the propagation direction of the field.

The problem can be considered as a 2D problem, due to symmetry. For this reason, I was thinking of implementing a 2D mesh and using a 3D vector electric field defined over the 2D mesh.

Before of implementing this kind of simulation, I tried to solve a simpler problem (the propagation of an electromagnetic plane wave) with this code:

from fenics import *

# Create mesh and define function space
mesh = UnitSquareMesh(12, 12)
V = FunctionSpace(mesh, 'N1curl', 2)

# Define k-vector
k_x = Constant(2*pi)
k_y = Constant(0)
k_z = Constant(0)
k = Constant((k_x, k_y, k_z))

# Define phase
phi = 0

# Define Field Polarization
u_x = 0
u_y = 1
u_z = 1

# Define the vector n normal on the boundary
n = FacetNormal(mesh)

# Define where to apply boundary conditions 
def boundary(x, on_boundary):
    return on_boundary

# Define the expression for the field on the boundary
u_D = Expression(('u_x*cos(k_x*x[0]+k_y*x[1]+phi)',
                'u_y*cos(k_x*x[0]+k_y*x[1]+phi)',
                'u_z*cos(k_x*x[0]+k_y*x[1]+phi)'),
                 degree=2, u_x=u_x, u_y=u_y, u_z=u_z, 
                 k_x=k_x, k_y=k_y, phi=phi, domain = mesh)

# Define the Dirichlet boundary condition
bc = DirichletBC(V, u_D, boundary)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
F = dot(k, k)*inner(u, v)*dx - inner(curl(u), curl(v))*dx
a, L = lhs(F), rhs(F)

# Compute solution
u = Function(V)
solve(a == L, u, bc)

That gives me this error:

keyTraceback (most recent call last):
  File "plane_wave.py", line 36, in <module>
    bc = DirichletBC(V, u_D, boundary)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/dirichletbc.py", line 131, in __init__
    super().__init__(*args)
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     fenics-support@googlegroups.com
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error:   Unable to create Dirichlet boundary condition.
*** Reason:  Illegal value dimension (3), expecting (2).
*** Where:   This error was encountered inside DirichletBC.cpp.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.dev0
*** Git changeset:  unknown
*** -------------------------------------------------------------------------

This error maybe can be solved easily, but it let me think about the feasibility of what I want to do: solving for a 3D vector field over a 2D geometry in fenics, especially when N1curl elements are involved.

Is this possible? Or are there any issues in doing that?

Thanks.

2

One should carefully consider the meaning of the H(\mathrm{curl};\Omega) conforming Nédélec element in the 2D and 3D setting. If you wish include some kind of slab symmetry, I expect the infinite axial component would belong to another space. You would then need to redefine the div, grad and curl UFL operators to accommodate this formulation.

I believe these techniques are investigated in Jian-Ming Jin’s book in the context of waveguides. For analysis I’d consult Peter Monk’s book. I vaguely recall it’s common practice when discretising waveguide models.

1 Like