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

**URL:** <https://fenicsproject.discourse.group/t/solving-for-a-3d-vector-field-over-a-2d-mesh-with-n1curl-elements/5216>\
**Category:** variational formulation\
**Created:** [March 3, 2021, 11:42am UTC](https://fenicsproject.discourse.group/t/solving-for-a-3d-vector-field-over-a-2d-mesh-with-n1curl-elements/5216 "2021-03-03T11:42:12Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![CastriMik](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/castrimik/32/2069_2.png) [@CastriMik](https://fenicsproject.discourse.group/u/CastriMik)\
**Post date:** [March 3, 2021, 11:42am UTC](https://fenicsproject.discourse.group/t/solving-for-a-3d-vector-field-over-a-2d-mesh-with-n1curl-elements/5216/1 "2021-03-03T11:42:12Z")

</div>

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:

```auto
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:

```auto
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.

---

<div class="post-metadata">

**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\
**Post date:** [March 3, 2021, 3:20pm UTC](https://fenicsproject.discourse.group/t/solving-for-a-3d-vector-field-over-a-2d-mesh-with-n1curl-elements/5216/2 "2021-03-03T15:20:10Z")

</div>

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](https://www.wiley.com/en-us/The+Finite+Element+Method+in+Electromagnetics%2C+3rd+Edition-p-9781118571361) in the context of waveguides. For analysis I’d consult Peter Monk’s [book](https://oxford.universitypressscholarship.com/view/10.1093/acprof:oso/9780198508885.001.0001/acprof-9780198508885). I vaguely recall it’s common practice when discretising waveguide models.