Hi all, I’m new to FEniCS, and as a first working example I wanted to simulate the reflection of a plane wave off a flat surface. I begin with a box-shaped Mesh with MeshTags indicating the reflective surface at y = 0, and the remaining boundaries. I define a TM plane wave function E_b coming in at a 45-degree angle. I want FEniCS to compute the ‘scattered’ function E_s, which I expect to be a plane wave travelling in the perpendicular direction.
With E = E_b + E_s, I use the double-curl equation
\qquad \nabla \times \nabla \times E - k^2 E = 0 \text{.}
At y = 0, I impose the boundary condition
\qquad n \times E = 0 \text{,}
while at the remaining boundaries, I impose the first-order absorbing boundary condition
\qquad n \times (\nabla \times E) + j k n \times (E \times n) = 0 \text{.}
I dump these equations straight into the code, as
a = ufl.inner(ufl.curl(ufl.curl(E_s)) - k**2 * E_s, v) * dx \
+ ufl.inner(ufl.cross(n, E_s), v) * ds_surface \
+ ufl.inner(ufl.cross(n, ufl.curl(E_s)) + 1j * k * ufl.cross(n, ufl.cross(E_s, n)), v) * ds_boundary
L = - ufl.inner(ufl.cross(n, E_b), v) * ds_surface
where
V = FunctionSpace(mesh, FiniteElement("N1curl", mesh.ufl_cell(), 3))
E_s = TrialFunction(V)
v = TestFunction(V)
dx = Measure("dx", my_mesh)
ds = Measure("ds", my_mesh, subdomain_data = facet_tags)
ds_surface = ds(1)
ds_boundary = ds(2)
n = FacetNormal(my_mesh)
I solve a == L as a LinearProblem, and the computation goes through without error, however the output of E_s is not remotely what I’m expecting.
Since I’m new to FEniCS (this is my first code snippet), I expect that I’m making a basic mistake in my setup. Does anyone have any suggestions?
Edit: Full code for completeness. I believe I’m using the dolfinx/lab:stable Docker image for running the code. The Git commit has for dolfinx is 24f86a9ce57df6978070dbee22b3eae8bb77235f.
import numpy as np
import pyvista as pv
import ufl
from dolfinx import fem, io, mesh, plot
from mpi4py import MPI
from petsc4py import PETSc
# Create box mesh
my_mesh = mesh.create_box(
MPI.COMM_WORLD,
[[0, 0, 0], [1, 1, 1]],
[10, 10, 1] # z-direction is irrelevant for now
)
boundaries = [(1, lambda x: np.isclose(x[1], 0)), # Reflective boundary
(2, lambda x: np.isclose(x[0], 1)) or np.isclose(x[0], 0) or np.isclose(x[1], 1) or np.isclose(x[2], 0) or np.isclose(x[2], 1)]
# Create tags for both boundaries
facet_indices, facet_markers = [], []
fdim = my_mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = mesh.locate_entities(my_mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tags = mesh.meshtags(my_mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
class BackgroundElectricField:
def __init__(self, theta: float, n_bkg: float, k0: complex):
self.theta = theta # Angle with respect to the x-axis!
self.n_bkg = n_bkg # Background refractive index
self.kx = self.n_bkg * k0 * np.cos(self.theta) # x wavevector
self.ky = self.n_bkg * k0 * np.sin(self.theta)
def eval(self, x):
phi = self.kx * x[0] + self.ky * x[1]
ax, ay = np.sin(self.theta), np.cos(self.theta)
return (-ax * np.exp(1j * phi), ay * np.exp(1j * phi), np.zeros(x[2].shape))
curl_el = ufl.FiniteElement("N1curl", my_mesh.ufl_cell(), 3) # ufl_cell or basix_cell depending on version
V = fem.FunctionSpace(my_mesh, curl_el)
theta = 1.75 * np.pi # Angle of incidence of the background field --- I think this is pointing to the down & right
n_bkg = 1.33 # Background refractive index
wavelength = 0.3 # Wavelength of the background field.
k = 2 * np.pi / wavelength
f = BackgroundElectricField(theta, n_bkg, k)
E_b = fem.Function(V)
E_b.interpolate(f.eval)
E_s = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
dx = ufl.Measure("dx", my_mesh)
ds = ufl.Measure("ds", my_mesh, subdomain_data = facet_tags)
ds_surface = ds(1)
ds_boundary = ds(2)
n = ufl.FacetNormal(my_mesh)
# Maxwell equations
a = ufl.inner(ufl.curl(ufl.curl(E_s)) - k**2 * E_s, v) * dx \
+ ufl.inner(ufl.cross(n, E_s), v) * ds_surface \
+ ufl.inner(ufl.cross(n, ufl.curl(E_s)) + 1j * k * ufl.cross(n, ufl.cross(E_s, n)), v) * ds_boundary
L = ufl.inner(ufl.cross(n, E_b), v) * ds_surface
# Solve problem
problem = fem.petsc.LinearProblem(a, L, bcs=[], petsc_options={"ksp_type": "preonly", "pc_type": "lu"})
E_s = problem.solve()
# Plot results
V_dg = fem.VectorFunctionSpace(my_mesh, ("Discontinuous Lagrange", 3))
E_dg = fem.Function(V_dg)
E_dg.interpolate(E_s)
pv.start_xvfb()
plotter = pv.Plotter()
plotter.set_position([0, 0, 0])
topology, cell_types, geometry = plot.create_vtk_mesh(V_dg)
grid = pv.UnstructuredGrid(topology, cell_types, geometry)
E_values = np.zeros((geometry.shape[0], 3), dtype = np.float64)
E_values[:, : my_mesh.topology.dim] = E_dg.x.array.reshape(geometry.shape[0], my_mesh.topology.dim).real # Discard imaginary parts for plotting
grid.point_data["E_sh"] = E_values
plotter = pv.Plotter()
plotter.add_mesh(grid.copy(), show_edges=True)
plotter.view_xy()
plotter.link_views()
plotter.show(jupyter_backend = 'panel')