Hello,
I am new to BEMPP and want to use it to simulate the scattering field of a dielectric cube backed by a PEC ground. The structure is illustrated in the figure below.
Following the two tutorial examples (flat screens and dielectric spheres), and the paper, I derived the equations using Maxwell’s boundary operators.
For the dielectric boundary:
$$(A_d^++A_d^-)V_d^++(A_{dp}^++A_{dp}^{-})V_p^+=V_d^{inc},$$
and for the PEC boundary:
$${\rm H_{pd}^+}(\gamma_t E_{d}^+)+{\rm E_{pd}^+}(\gamma_N E_{d}^+)+{\rm E_{p}^+}(\gamma_N E_{p}^+)=\gamma_tE_p^{inc}$$
I followed the Laplace mixed-Neumann Dirichlet example to define the function spaces, boundary operators, and grid functions for the incident field. And I used the LU solver. However, the solutions seem incorrect. Some potential values around the cube are extremely large, as shown in the left figure below. As a comparison, I simulate the flat screen case without the dielectric cube, as shown in the right figure below. The result seems correct with specular reflection.
I have attached the code below. I feel that the equations I derived have problems. Then how should I assemble the boundary operators for this structure? I would appreciate any suggestions.
import bempp.api
from bempp.api.operators.boundary.maxwell import electric_field, magnetic_field
import numpy as np
import time
from matplotlib import pyplot as plt
from mpl_toolkits.axes_grid1 import make_axes_locatable
from matplotlib.patches import Rectangle
bempp.api.enable_console_logging()
bempp.api.DEFAULT_DEVICE_INTERFACE = "opencl"
bempp.api.DEFAULT_PRECISION = 'double'
bempp.api.PLOT_BACKEND = "paraview"
start = time.time()
obj_size = 4
obj_thickness = 1
grid = bempp.api.shapes.shapes.cuboid(length=(obj_size,obj_size,obj_thickness), origin=(-obj_size/2,-obj_size/2,0), h=0.1)
dielectric_segments = [2,3,4,5,6]
metal_segments = [1]
wavelength = 1
eps_0 = 8.854187817E-12
mu_0 = 4 * np.pi * 1E-7
eps_r = 2.6
mu_r = 1.0
wavelength_r = wavelength / np.sqrt(eps_r*mu_r)
k_ext = 2 * np.pi / wavelength
k_int = 2 * np.pi / wavelength_r
theta_inc = np.pi/4
phi_inc = 0
direction = np.array([np.cos(phi_inc)*np.sin(theta_inc), np.sin(phi_inc)*np.sin(theta_inc), np.cos(theta_inc)])
polarization = np.array([-np.sin(phi_inc), np.cos(phi_inc), 0])
def plane_wave(point):
return polarization * np.exp(-1j * k_ext * np.dot(point, direction))
@bempp.api.complex_callable
def tangential_trace(point, n, domain_index, result):
value = polarization * np.exp(-1j * k_ext * np.dot(point, direction))
result[:] = np.cross(value, n)
@bempp.api.complex_callable
def neumann_trace(point, n, domain_index, result):
value = -np.cross(direction, polarization) * 1j * k_ext * np.exp(-1j * k_ext * np.dot(point, direction))
result[:] = 1.0 / (1j * k_ext) * np.cross(value, n)
global_div_space = bempp.api.function_space(grid, "RWG", 0)
global_curl_space = bempp.api.function_space(grid, "SNC", 0)
div_space_dielectric_segment = bempp.api.function_space(grid, "RWG", 0, segments=dielectric_segments)
div_space_metal_segment = bempp.api.function_space(grid, "RWG", 0, segments=metal_segments)
curl_space_dielectric_segment = bempp.api.function_space(grid, "SNC", 0, segments=dielectric_segments)
curl_space_metal_segment = bempp.api.function_space(grid, "SNC", 0, segments=metal_segments)
H_DD_int = magnetic_field(div_space_dielectric_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_int,assembler="dense")
E_DD_int = electric_field(div_space_dielectric_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_int,assembler="dense")
H_DD_ext = magnetic_field(div_space_dielectric_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_ext,assembler="dense")
E_DD_ext = electric_field(div_space_dielectric_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_ext,assembler="dense")
H_DP_int = magnetic_field(div_space_metal_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_int,assembler="dense")
E_DP_int = electric_field(div_space_metal_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_int,assembler="dense")
H_DP_ext = magnetic_field(div_space_metal_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_ext,assembler="dense")
E_DP_ext = electric_field(div_space_metal_segment,div_space_dielectric_segment,curl_space_dielectric_segment,k_ext,assembler="dense")
H_PD_ext = magnetic_field(div_space_dielectric_segment,div_space_metal_segment,curl_space_metal_segment,k_ext,assembler="dense")
E_PD_ext = electric_field(div_space_dielectric_segment,div_space_metal_segment,curl_space_metal_segment,k_ext,assembler="dense")
E_PP_ext = electric_field(div_space_metal_segment,div_space_metal_segment,curl_space_metal_segment,k_ext,assembler="dense")
A = bempp.api.GeneralizedBlockedOperator([[H_DD_ext + H_DD_int, E_DD_ext + E_DD_int, E_DP_ext + E_DP_int], \
[-(E_DD_ext + E_DD_int), H_DD_ext + H_DD_int, H_DP_ext + H_DP_int], \
[H_PD_ext, E_PD_ext, E_PP_ext]])
E_DD_inc_t = bempp.api.GridFunction(div_space_dielectric_segment, fun=tangential_trace, dual_space=curl_space_dielectric_segment)
E_DD_inc_n = bempp.api.GridFunction(div_space_dielectric_segment, fun=neumann_trace, dual_space=curl_space_dielectric_segment)
E_PP_inc_t = bempp.api.GridFunction(div_space_metal_segment, fun=tangential_trace, dual_space=curl_space_metal_segment)
rhs = [E_DD_inc_t, E_DD_inc_n, E_PP_inc_t]
(sol_DD_ext_t, sol_DD_ext_n, sol_PP_ext_n) = bempp.api.linalg.lu(A, rhs)
end = time.time()
print(end-start)
# plot near field
nx = 200
nz = 200
x, y, z = np.mgrid[-5:5:nx*1j, 0:0:1j, -5:5:nz*1j]
points = np.vstack((x.ravel(), y.ravel(), z.ravel()))
interior_indices = np.logical_and.reduce((points[0, :] > -obj_size/2, points[0, :] < obj_size/2, \
points[1, :] > -obj_size/2, points[1, :] < obj_size/2, \
points[2, :] > 0, points[2, :] < obj_thickness))
exterior_indices = ~interior_indices
ext_points = points[:, exterior_indices]
int_points = points[:, interior_indices]
hpot_DD_int = bempp.api.operators.potential.maxwell.magnetic_field(sol_DD_ext_t.space, int_points, k_int, assembler="dense")
epot_PP_int = bempp.api.operators.potential.maxwell.electric_field(sol_PP_ext_n.space, int_points, k_int, assembler="dense")
epot_DD_int = bempp.api.operators.potential.maxwell.electric_field(sol_DD_ext_n.space, int_points, k_int, assembler="dense")
hpot_DD_ext = bempp.api.operators.potential.maxwell.magnetic_field(sol_DD_ext_t.space, ext_points, k_ext, assembler="dense")
epot_PP_ext = bempp.api.operators.potential.maxwell.electric_field(sol_PP_ext_n.space, ext_points, k_ext, assembler="dense")
epot_DD_ext = bempp.api.operators.potential.maxwell.electric_field(sol_DD_ext_n.space, ext_points, k_ext, assembler="dense")
exterior_values = -hpot_DD_ext * sol_DD_ext_t - epot_DD_ext * sol_DD_ext_n - epot_PP_ext * sol_PP_ext_n
interior_values = hpot_DD_int * sol_DD_ext_t + epot_DD_int * sol_DD_ext_n + epot_PP_int * sol_PP_ext_n
plt.rcParams['figure.figsize'] = (20, 16)
scattered_field = np.empty((3, points.shape[1]), dtype='complex128')
scattered_field[:, :] = np.nan
scattered_field[:, exterior_indices] = exterior_values
scattered_field[:, interior_indices] = interior_values
total_field = np.empty((3, points.shape[1]), dtype='complex128')
for ext_ind in np.arange(points.shape[1])[exterior_indices]:
total_field[:, ext_ind] = scattered_field[:, ext_ind] + plane_wave(points[:, ext_ind])
total_field[:, interior_indices] = interior_values
squared_scattered_field = np.sum(np.abs(scattered_field)**2, axis=0)
squared_total_field = np.sum(np.abs(total_field)**2, axis=0)
scattered_image = squared_scattered_field.reshape(nx, nz).T
total_image = squared_total_field.reshape(nx, nz).T
fig, axes = plt.subplots(1, 2, sharex=True, sharey=True)
f0 = axes[0].imshow(scattered_image, origin='lower', cmap='magma', extent=[-4, 4, -4, 4], vmin=0, vmax=1000)
axes[0].add_patch(
Rectangle((-obj_size/2, 0), obj_size, obj_thickness, facecolor='None', edgecolor='k')
)
axes[0].set_title("Squared scattered field strength")
divider = make_axes_locatable(axes[0])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(f0, cax=cax)
f1 = axes[1].imshow(total_image, origin='lower', cmap='magma', extent=[-4, 4, -4, 4], vmin=0, vmax=1000)
axes[1].add_patch(
Rectangle((-obj_size/2, 0), obj_size, obj_thickness, facecolor='None', edgecolor='k')
)
axes[1].set_title("Squared total field strength")
divider = make_axes_locatable(axes[1])
cax = divider.append_axes("right", size="5%", pad=0.05)
fig.colorbar(f1, cax=cax)
plt.show()
# plot far field
number_of_angles = 400
angles = np.pi * np.linspace(0, 1, number_of_angles)
unit_points = np.array([np.cos(angles), np.zeros(number_of_angles), np.sin(angles)])
far_field = np.zeros((3, number_of_angles), dtype='complex128')
hfar_DD_ext = bempp.api.operators.far_field.maxwell.magnetic_field(sol_DD_ext_t.space, unit_points, k_ext, assembler="dense")
efar_PP_ext = bempp.api.operators.far_field.maxwell.electric_field(sol_PP_ext_n.space, unit_points, k_ext, assembler="dense")
efar_DD_ext = bempp.api.operators.far_field.maxwell.electric_field(sol_DD_ext_n.space, unit_points, k_ext, assembler="dense")
far_field = -hfar_DD_ext * sol_DD_ext_t - efar_DD_ext * sol_DD_ext_n - efar_PP_ext * sol_PP_ext_n
plt.rcParams['figure.figsize'] = (10, 8)
cross_section = 10 * np.log10(4 * np.pi * np.sum(np.abs(far_field)**2, axis=0))
plt.plot(angles * 180 / np.pi, cross_section)
plt.title("Scattering Cross Section [dB]")
_ = plt.xlabel('Angle (Degrees)')
plt.show()