### Maxwell Equations Boundary Conditions

291
views
1
6 months ago by

Hello,

since my previous post was a little to stumpy I will retry to describe my problem.

I want to solve Maxwells equations for a scattering problem. The underlying geometry is a small metallic (perfect conductor) square with surface $S_2$S2 which is surrounded by a larger square with surface $S_1$S1 which acts as an absorbing boundary. The inner of the larger square minus the smaller square is the area  $\Omega$Ω . I want to calculate the scattered field  $\vec{E}^{sc}$Esc knowing the incoming field  $\vec{E}^{in}$Ein . The functional for this problem is given by:  $F\left(\vec{E}^{sc}\right)=F_1\left(\vec{E}^{sc}\right)+F_2\left(\vec{E}^{sc}\right)$F(Esc)=F1(Esc)+F2(Esc) with $F_1\left(\vec{E}^{sc}\right)=\int_{\Omega}\nabla\times\vec{E}^{sc}\cdot\nabla\times\vec{E}^{sc}dA+k_0^2\int_{\Omega}\vec{E}^{sc}\cdot\vec{E}^{sc}dA$F1(Esc)=Ω×Esc·×EscdA+k02ΩEsc·EscdA and   $F_2\left(\vec{E}^{sc}\right)=i\cdot k_0\cdot\int_{S_1}\vec{n}\times\vec{E}^{sc}\cdot\vec{n}\times\vec{E}^{sc}ds$F2(Esc)=i·k0·S1n×Esc·n×Escds . The first part accounts for the integration of the Helmholtz equation and the second parts accounts for the absorbing boundary condition. The solution can now be found by minimizing  $F\left(\vec{E}^{sc}\right)$F(Esc) under the condition that $\vec{n}\times\vec{E}^{sc}=-\vec{n}\times\vec{E}^{in}$n×Esc=n×Ein at $S_2$S2. My problem is to implement the boundary condition at $S_2$S2.

Any suggestions or ideas are welcome :) Thanks for your help!

from fenics import *
from mshr import *
import numpy as np
import scipy.sparse as sps
import matplotlib.pyplot as plt

from linear_algebra import convert_matrix_dolfin_scipy, convert_vector_dolfin_numpy, solve_scipy_system

"""
Define natural constants
"""
c0 = 299792458
mu0 = 4*np.pi*1e-7
eps0 = 1/(c0**2*mu0)

"""
Define Problem parameters
"""
# physical properties
eps_r = 1
mu_r = 1
freq = 1e9
lam_0 = c0/freq
k_0 = 2*np.pi*lam_0
E_0 = 1.
r_0 = np.array([0, -10.])
k = k_0 * np.array([1., 0.])

"""
Define parameters for solvers and backends
"""
parameters['linear_algebra_backend'] = 'Eigen'

"""
Import the needed geometry
"""
mesh = Mesh('geometry/square/square.xml')
domains = MeshFunction('size_t', mesh, 'geometry/square/square_physical_region.xml')
boundaries = MeshFunction('size_t', mesh, 'geometry/square/square_facet_region.xml')

"""
Define incoming wave as matrix expression
E = E_0*p*exp(-i*k*(r-r_0)), p = cross(e_r, e_z)
with (((Re(E_x), Re(E_y)), ((Im(E_x), Im(E_y)))
"""
E_in = Expression(((('E_0 * (x[1] - r_0_y) * cos(- k_x * (x[0] - r_0_x) - k_y * (x[1] - r_0_y)) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))'),\
('E_0 * (x[1] - r_0_y) * sin(- k_x * (x[0] - r_0_x) - k_y * (x[1] - r_0_y)) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))')),\
(('-E_0 * (x[0] - r_0_x) * cos(- k_x * (x[0] - r_0_x) - k_y * (x[1] - r_0_y)) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))'),\
('-E_0 * (x[0] - r_0_x) * sin(- k_x * (x[0] - r_0_x) - k_y * (x[1] - r_0_y)) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))'))),\
E_0 = E_0, k_x = k[0], k_y = k[1], r_0_x = r_0[0], r_0_y = r_0[1], degree = 1, domain = mesh.ufl_domain())
E_in_abs = Expression((('E_0 * (x[1] - r_0_y) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))'),\
('-E_0 * (x[0] - r_0_x) / sqrt( pow(x[0] - r_0_x, 2) + pow(x[1] - r_0_y, 2))')),\
E_0 = E_0, r_0_x = r_0[0], r_0_y = r_0[1], degree = 1, domain = mesh.ufl_domain())

"""
Define function space for test and trial function
"""
# Function space
func_space = FunctionSpace(mesh, 'N1curl', 1)
# Trial and Test function
trial_func = TrialFunction(func_space)
test_func = TestFunction(func_space)
# Measures
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Normal vector to facet
face_norm = FacetNormal(mesh)

"""
Define the bilinear forms, assemble and convert to scipy matrix
"""
bilform_curl = dot(curl(trial_func), curl(test_func))*dx
bilform_dot = dot(trial_func, test_func)*dx
bilform_face = (trial_func[0]*face_norm[1]-trial_func[1]*face_norm[0])*(test_func[0]*face_norm[1]-test_func[1]*face_norm[0])*ds(1)

# Assemble the matrices
mat_curl = assemble(bilform_curl)
mat_dot = assemble(bilform_dot)
mat_face = assemble(bilform_face)

# Convert to scipy sparse
mat_curl = convert_matrix_dolfin_scipy(mat_curl)
mat_dot = convert_matrix_dolfin_scipy(mat_dot)
mat_face = convert_matrix_dolfin_scipy(mat_face)

mat = 1/mu_r*mat_curl - k_0**2*eps_r*mat_dot + 1j*k_0*mat_face

Community: FEniCS Project

Hello,

Am I asking the question wrong? Is there any information missing? I really need some help to proceed with the simulation.

Thanks - Christopher

written 6 months ago by Christopher Zachow
1
DirichletBC class for boundary conditions fixes values of degrees of freedom. Just see what DOFs of N1curl space are (in periodic table of finite elements or directly in the fenics source code :( ) and you'll find out, that by specifying the dirichletbc on a boundary you are effectively fixing $n \times E$.
written 6 months ago by Michal Habera

Since I am solving a 2D problem there is no natural way to express the cross product. If I use the 'DirichletBC' it expects a 2-vector not a value for the z-component.

written 6 months ago by Christopher Zachow

2
6 months ago by
You must supply vector $E$ on the boundary to the dirichlet BC. It will do 2D cross product $n \times E$, which is a number for vectors in 2D and this number is used as the dof-value.

Mathematical note: cross (exterior) product in 2D makes sense, but the result doesn't stay in the same exterior algebra (it is 0-form, not 1-form)...

Play with DirichletBC's method get_boundary_values() (it returns global degree of freedom index and value of that dof) and you'll see it does what I said.
Finally it works, thanks a lot!
written 6 months ago by Christopher Zachow